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Abstract Oxygen is a promising exoplanet biosignature due to the evolutionary advantage conferred 
by harnessing starlight for photosynthesis, and the apparent low likelihood of maintaining oxygen- 

rich atmospheres without life. Hypothetical scenarios have been proposed for non-biological oxygen 
accumulation on planets around late M-dwarfs, where the extended pre-main sequence may favor abiotic 
O, accumulation. In contrast, abiotic oxygen accumulation on planets around F, G, and K-type stars is 
seemingly less likely, provided they possess substantial non-condensable gas inventories. The comparative 
robustness of oxygen biosignatures around larger stars has motivated plans for next-generation 
telescopes capable of oxygen detection on planets around sun-like stars. However, the general tendency 
of terrestrial planets to develop oxygen-rich atmospheres across a broad range of initial conditions and 
evolutionary scenarios has not been explored. Here, we use a coupled thermal-geochemical-climate 
model of terrestrial planet evolution to illustrate three scenarios whereby significant abiotic oxygen can 
accumulate around sun-like stars, even when significant noncondensable gas inventories are present. For 
Earth-mass planets, we find abiotic oxygen can accumulate to modern levels if (1) the CO;:H5O ratio of 
the initial volatile inventory is high, (2) the initial water inventory exceeds ~50 Earth oceans, or (3) the 
initial water inventory is very low («0.3 Earth oceans). Fortunately, these three abiotic oxygen scenarios 
could be distinguished from biological oxygen with observations of other atmospheric constituents 

or characterizing the planetary surface. This highlights the need for broadly capable next-generation 
telescopes that are equipped to constrain surface water inventories via time-resolved photometry and 
search for temporal biosignatures or disequilibrium combination biosignatures to assess whether oxygen 
is biogenic. 


Plain Language Summary Next-generation telescopes will search for life on exoplanets by 
looking for the spectral signatures of biogenic gases. Oxygen has been considered a reliable biosignature 
gas, especially for planets around sun-like stars where non-biological, photochemical production is 
unlikely. This motivates plans for future telescopes specifically designed for oxygen detection. Here, we 
develop a coupled model of the atmosphere-interior evolution of terrestrial planets to show that lifeless 
planets in the habitable zone could develop oxygen-rich atmospheres relatively easily. These false positives 
for biological oxygen could be distinguished from inhabited planets using other contextual clues, but their 
existence implies next-generation telescopes need to be capable of characterizing planetary environments 
and searching for multiple lines of evidence for life, not merely oxygen. 


1. Introduction 


The search for life beyond Earth is a key motivator for exoplanet astronomy. Among the various life detec- 
tion approaches that have been proposed, atmospheric oxygen is arguably the most promising biosignature. 
This is because any organism that adapts to exploit free energy from starlight will have a competitive advan- 
tage over organisms that are limited by geochemical sources of free energy. The earliest incontrovertible ev- 
idence for life on Earth around 3.5 Ga coincides with fossilized photosynthetic stromatolites (Buick, 2008). 
While it is unknown whether the evolution of oxygenic photosynthesis is contingent—the metabolism 
emerged only once in Earth's history (Mulkidjanian et al., 2006)—oxygenic photosynthesis confers a unique 
evolutionary advantage over other forms of photosynthesis since the required substrates, carbon dioxide 
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(CO;) and water (H5O) are likely ubiquitous on habitable worlds. Moreover, the accumulation of biogenic 
oxygen in planetary atmospheres, as has occurred on Earth (Holland, 2006; Lyons et al., 2014), is readily 
detectable over interstellar distances thanks to absorption features in the visible, near IR, and (for ozone) 
thermal infrared (reviewed by Meadows et al., 2018). Indeed, three of the four proposed mission concepts 
currently under consideration by the 2020 Astrophysics Decadal Survey are specifically designed to be ca- 
pable of life detection via oxygen or ozone (Fischer et al., 2019; Gaudi et al., 2020; Meixner et al., 2019). 
Ground-based Extremely Large Telescopes may also be capable of oxygen detection via high resolution 
spectroscopy (Leung et al., 2020; Rodler & López-Morales, 2014; Snellen et al., 2013). The James Webb 
Space Telescope (JWST) could conceivably detect oxygen/ozone biosignatures for nearby transiting planets, 
although this would probably require a prohibitively large number of transits for known targets (Fauchez 
et al., 2020; Krissansen-Totton, Garland, et al., 2018; Lustig-Yaeger et al., 2019; Wunderlich et al., 2019). 
Finally, oxygen is often considered to be a good biosignature because—apart from the exceptions discussed 
below—it is seemingly difficult for habitable zone terrestrial planets to maintain oxygen-rich atmospheres 
without life (reviewed by Meadows et al., 2018). 


The terrestrial planets that will be accessible to spectroscopic characterization with JWST and ground-based 
ELTs will orbit M-dwarfs due to the favorable signal-to-noise of M-dwarf transits. This is unfortunate from 
the standpoint of recognizing oxygen biosignatures because several features of M-dwarfs make them suscep- 
tible to non-biological oxygen accumulation. In particular, the extended pre-main sequence of late M-dwarfs 
could yield habitable zone terrestrial planets with hundreds or thousands of bar O; from XUV-driven hy- 
drogen loss (Luger & Barnes, 2015). At least some of this oxygen will likely dissolve in a surface magma 
ocean and be sequestered in the mantle, but retaining oxygen-rich atmospheres is still possible, especially 
for highly irradiated terrestrial planets (Barth et al., 2020; Schaefer et al., 2016; Wordsworth et al., 2018). The 
early development of abiotic, oxygen-rich atmospheres could prevent the subsequent emergence of life by 
precluding prebiotic chemistry (Wordsworth et al., 2018). It has also been argued that the spectral energy 
distributions of M-dwarfs are favorable for the photochemical accumulation of O, and O; (Domagal-Gold- 
man et al., 2014; Gao et al., 2015; Harman et al., 2015; Tian et al., 2014). While recent assessments of the 
role of lightning in such photochemical models (Harman et al., 2018) and improved near-UV water cross 
sections (Ranjan et al., 2020) may preclude some of these scenarios, photochemical runaways yielding O3- 
CO rich atmospheres remain a strong possibility for late M dwarfs (Hu et al., 2020). More fundamentally, 
however, whether M-dwarf terrestrial planets can sustain habitable surface conditions, or indeed any kind 
of substantial atmosphere, for billions of years given the harsh stellar radiation environment is unknown 
(Airapetian et al., 2017; Dong et al., 2018; Shields et al., 2016; Zahnle & Catling, 2017). Premature claims of 
biogenic oxygen on M-dwarf planets will rightfully be met with skepticism. 


These challenges to M-dwarf habitability and potential ambiguities with oxygen biosignatures motivates 
the need to characterize habitable zone planets around F, G, and K-type stars (Arney et al., 2019; National 
Academies of Sciences & Medicine, 2018). The relatively short pre-main sequence of larger stars mean 
that the accumulation of non-biological oxygen is less likely. Indeed, the most plausible mechanism for 
non-biological O, accumulation on habitable zone planets around sun-like stars is due to small atmospher- 
ic inventories of non-condensable species like N»: this weakens the water cold trap and may produce high 
stratospheric water abundances and therefore elevated H escape rates (Kleinbóhl et al., 2018; Wordsworth & 
Pierrehumbert, 2014). Whether terrestrial planets are likely to form with low nitrogen atmospheric invento- 
ries is unknown, but in principle, it might be possible to rule out low non-condensable scenarios by looking 
for the spectral signatures of abundant nitrogen (Schwieterman et al., 2015) or via inferring background gas 
abundances from retrieved surface pressure (Feng et al., 2018). 


The broader question as to whether habitable zone terrestrial planet evolution tends toward anoxic atmos- 
pheres, or whether non-biological oxygen production via water photodissociation and hydrogen escape can 
overwhelm geochemical sinks, has not been explored. The robustness of oxygen biosignatures rests on the 
assumption that for temperate planets with effective cold traps, small abiotic oxygen source fluxes from H 
escape will be overwhelmed by geological sinks. To test this assumption, it is necessary to model the redox 
evolution of terrestrial planet from formation onwards. This is because planetary redox evolution depends 
on both the initial state of the atmosphere and mantle after the magma ocean has solidified, and on the 
subsequent internal evolution and atmospheric state. Interior evolution dictates crustal production rates 
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Figure 1. Schematic of terrestrial planet geochemical evolution model. The planetary redox budget, thermal-climate 
evolution, and volatile budget are modeled from initial magma ocean (left) through to temperate geochemical cycling 
(right). Oxygen fluxes are shown by green arrows, energy fluxes by black arrows, carbon fluxes by orange arrows, 

and water fluxes by blue arrows. Note that the net loss of H to space adds oxygen to the atmosphere. During the 
magma ocean phase, the radius of solidification, rs, begins at the core-mantle boundary and moves toward the surface 
as internal heat is dissipated. The rate at which this occurs is controlled by radiogenic heat production, Qradioactives 
convective heatflow from the mantle to the surface, dmantie, and core heatflow, Qore. This internal heatflow balances the 
difference between outgoing longwave radiation, OLR, and incoming shortwave radiation, ASR. The oxygen fugacity of 
the mantle, fO;, and the water and carbon content mantle and surface reservoirs are tracked throughout. 
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and outgassing fluxes, which determine the efficiency of geologic sinks of oxygen. Moreover, self-consistent 
modeling of surface climate and weathering processes is necessary because these modulate surface volatile 


inventories and oxygen production via atmospheric H escape. 


Here, we explore the tendency of terrestrial planets to produce abiotic, oxygen-rich atmospheres with a 
fully coupled model of planetary redox, climate, and thermal evolution. Given the vagaries of planetary 
accretion, migration, and core formation, it is likely that habitable zone terrestrial planets exist with a broad 
range of initial volatile inventories (Raymond et al., 2004, 2013; Righter, 2015). Indeed, both planet forma- 
tion models (Raymond et al., 2004) and exoplanet demographics suggest that water-rich terrestrial planets 
are abundant (Zeng et al., 2019), and several candidate waterworlds are amenable to JWST characterization 
(Benneke et al., 2019; Grimm et al., 2018). In contrast, transit and thermal phase curve observations show 
that the terrestrial exoplanet LHS 3844b is either a bare rock or possesses a thin, high mean-molecular 
weight atmosphere susceptible to stellar wind erosion (Diamond-Lowe et al., 2020; Kreidberg et al., 2019), 
and modeling implies it must have formed with less volatiles than the Earth (Kane et al., 2020). Relative 
abundances of volatiles may also be variable; C/H mass ratios in carbonaceous chondrites vary from ~1.6 to 
6.4 (Nittler et al., 2004). We thus apply our model to explore the possibility of non-biological oxygen accu- 
mulation for a range of initial volatile inventories. Crucially, we avoid varying background N, abundances 
because our focus is on investigating oxygen false positives on planets with comparatively high background 
non-condensables inventories, that is, sufficient to maintain a cold trap under temperate conditions. 


2. Materials and Methods 


Our model is summarized schematically in Figure 1. A complete description is provided in the supporting 
information and the Python code is open source (https://doi.org/10.5281/zenodo.4539040). Here, we sum- 
marize the salient features of the model necessary for interpreting key findings. There are no biological 
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sources and sinks of volatiles because we are testing the capacity of lifeless planets to accumulate oxygen. 
Unless stated otherwise, all calculations assume an Earth-mass planet at 1 AU. 


Planetary evolution is divided into an initial magma ocean phase, and a subsequent solid-mantle phase, 
as shown in Figure 1, although a planet may transition between magma ocean and solid mantle multiple 
times. The model is initialized with a fully molten mantle and some endowment of volatiles, radiogenic 
inventory, and an initial mantle oxygen fugacity (i.e., after core formation, Figure 1, left). Magma ocean 
solidification follows previous models such as Lebrun et al. (2013) and Schaefer et al. (2016): the magma 
ocean freezes from the core, upwards, as governed by the following equation: 


2 2 dr, dT, 
Vinantle Pn O radioactive x 474, * Deore + 470,5, H fusions FS E VinantePm€ p dt (1) 


Here, Vmante is the volume of the molten mantle, p, is the average density of the mantle, Qradioactive is radio- 
genic heat production per unit mass, r, is planetary radius, Hfusion is the latent heat of fusion of silicates, rs 
is the solidification radius, c; is the specific heat of silicates, Qeore is the heatflow from the metallic core, and 
T, is mantle potential temperature. The heatflow from the interior, qm, is calculated using 1-D convective 
parameterization, with temperature-dependent magma ocean viscosity, v(T;), (see Figure S3 and support- 
ing information Section A.3): 


(2) 


Here, Tsurp is mean surface temperature, and C, is a constant that depends on thermal conductivity, thermal 
diffusivity, critical Rayleigh number, gravity, and thermal expansivity. Supporting information Section A.1 
describes this convection parameterization in more detail. Equation 1 continues to govern the thermal evo- 
lution of the mantle after the magma ocean has solidified with dr,/dt = 0.0. 


During magma ocean solidification H, C, and O are partitioned between dissolved melt phases, crystalline 
phases, and the atmosphere by assuming chemical equilibrium (supporting information Section A.7). The 
rate at which the mantle freezes is controlled by outgoing longwave radiative (OLR), which is balanced by 
heat from interior and absorbed shortwave radiation (ASR) from the host star at every timestep (see climate 
model description below): 


dm + ASR = OLR (3) 


During the magma ocean phase, planetary oxidation may occur from the loss of hydrogen to space (less 
oxygen drag): 


H,0 > H; (space) + 0.50, (4) 


Free oxygen produced via H escape is dissolved in the melt and may be transferred to the solid mantle as 
the magma ocean solidifies (Schaefer et al., 2016). We parameterized atmospheric escape as either diffu- 
sion limited or XUV-limited, depending on the composition of the stratosphere and the stellar XUV flux 
(supporting information Section A.9). During XUV-driven escape of a steam-dominated atmosphere, the 
hydrodynamic escape of H may drag along O (and even CO,) following Odert et al. (2018). In contrast, if 
the stratosphere is mostly dry, then the escape of H will be limited by the rate at which H,O can diffuse 
through the cold trap (Wordsworth & Pierrehumbert, 2013), and nothing heavier than H escapes. Standard 
parameterizations of solar bolometric luminosity (Baraffe et al., 1998, 2002) and XUV luminosity evolution 
are adopted (Tu et al., 2015), as described in supporting information Section A.4. 


A radiative-convective climate model is used to self-consistently calculate surface temperature, OLR, ASR, 
the water vapor profile, and surface liquid water inventory (if any) during both the magma ocean phase 
and subsequent temperate evolution. OLR is a function of the surface H,O and CO, inventories, and is 
calculated using the publicly available correlated-k radiative transfer code of Marcq et al. (2017). To obtain 
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OLR in the presence of condensable water vapor, a dry adiabat to moist adiabat to isothermal atmospheric 
structure is assumed (Kasting, 1988). To calculate ASR across a wide range of temperatures, we adapted the 
albedo parameterization described in Pluriel et al. (2019). Refer to supporting information Section A.5 for 
full details of radiative transfer calculations along with example outputs. 


When heat from accretion and short-lived radiogenics is sufficiently dissipated—the timescale for which is 
controlled by insolation and greenhouse warming from outgassed volatiles—a planet's surface temperature 
may drop below the solidus and the magma ocean phase is over. At this point, the model transitions to sol- 
id-state mantle convection and temperate geochemical cycling (Figure 1, right). The redox budget during 
solid-state evolution is modeled as follows: the only source of oxygen is still atmospheric escape using the 
same parameterization as described above. However, there are now numerous crustal sinks for oxygen in- 
cluding (i) subaerial and submarine outgassing of reduced species (e.g., H}, CO, and CH,), (ii) water-rock 
serpentinizing reactions that generate H, (the “wet crustal” sink), and (iii) direct oxidation of surface crust 
by atmospheric oxygen (the “dry crustal” sink): 


H, + 0.50, > H,O 
üh CO +0.50, > CO, 
CH, +20, > CO, + 2H,0 (5) 
(ii),2FeO + H,O > H, + 2FeO, ; 
(iii), 2FeO + 0.50, — H, + 2FeO, ; 


The sizes of these three oxygen sinks are self-consistently calculated from the planetary interior evolution 
and mantle volatile content. Outgassing fluxes are calculated using the melt-gas equilibrium outgassing 
model of Wogan et al. (2020); outgassing fluxes depend on mantle oxygen fugacity, degassing overburden 
pressure, the volatile content of the mantle, specifically H;O and CO, content, and the rate at which melt 
(new crust) is produced (see supporting information Section A.10). The possible influence of pressure over- 
burden on redox evolution has been discussed previously (Wordsworth et al., 2018), and is quantifiable 
within our outgassing model framework. Dry and wet crustal sinks for oxygen similarly depend on crustal 
production rates, and are described in full in supporting information Sections A.12 and A.13. Crustal pro- 
duction is calculated from interior heatflow, which is modulated by temperature-dependent mantle viscos- 
ity and radiogenic heat production (supporting information Section A.10). We assume plate tectonics when 
calculating crustal production rates to maximize crustal sinks of oxygen. 


Weathering processes (Krissansen-Totton, Arney, et al., 2018) and the deep hydrological cycle (Schaefer & 
Sasselov, 2015) are explicitly modeled because climate and surface volatile inventories control crustal ox- 
ygen sinks and atmospheric escape processes. Our model incorporates a rudimentary carbon cycle, which 
is a simplified version of that described in Krissansen-Totton, Arney, et al. (2018). Carbon is added to the 
atmosphere via magmatic outgassing (described above), and returned via continental weathering and sea- 
floor weathering, whose relative contributions depend on climate and the total surface water inventory. Our 
carbon cycle parameterization is described in supporting information Section A.11 and the deep hydrologi- 
cal cycle in supporting information Section A.12. Self-consistent climate modeling enables an assessment of 
whether abiotic oxygen can coexist with a habitable surface climate, which would be especially problematic 
for unambiguous biosignature gas interpretations. 


The time evolution of volatile reservoirs during both the magma ocean phase and the solid state evolution 
is governed by the following equations (see supporting information Section A.6 for details): 


AM soia. di 
DAD. 2 as 
-- Amp nk», fra; Ts + Fingas—A; E Fougas-A; 
dt dt 
dM Fyid—A dr, (6) 
uid-A; — 2. s. 
dt = 4m nk, Sin dt ~ Fingas—A; + Foutgas- Aj ~ Esc, 


Here, A; represents a generic volatile species (e.g., H;O, CO,, free O). The first term on the right hand side 
represents the transfer of volatiles from the fluid phases (magma ocean + atmosphere) to the solid mantle 
as the magma ocean solidifies: k,, is the melt-solid partition coefficient for species A;, and fry, represents 
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the mass fraction of the volatile species in the magma. The remaining fluxes are subaerial plus submarine 
outgassing from the mantle to the atmosphere, Fouigas_a,, ingassing from the atmosphere to the mantle (e.g., 
crustal oxidation or hydration), Fingas-;, and escape to space Esc,,. Because we are assuming a plate tec- 
tonics regime, the model does not separately track volatile reservoirs in the crust and mantle. Instead, we 
assume carbon, oxygen, and water added to crust is immediately subducted into the mantle; a single, well- 
mixed “interior” reservoir is used to represent storage of volatiles in solid silicates. The extremely broad 
range of crustal hydration and crustal oxidation efficiency factors sampled (see below) can accommodate 
differing subduction and arc volcanism efficiencies. 


Note that the model only tracks C, H, and O-bearing species, as well as Fe**/Fe** speciation in the interior. 
Nitrogen fluxes are not modeled, and we instead assumed a 1 bar N, background partial pressure in all mod- 
el runs. This conservative assumption ensures that, for temperature surface conditions, there are always suf- 
ficient non-condensables to maintain a cold trap and prevent excessive water loss; any oxygen accumulation 
that results is due to other processes. 


The model does not include any explicit photochemistry; it tracks fluxes of oxygen into/out of the combined 
atmosphere-ocean reservoir, and all outgassed reductants are assumed to instantaneously deplete atmos- 
pheric oxygen. This simplification is adequate for estimating oxygen accumulation because if oxygen sourc- 
es exceed oxygen sinks, then oxidant build-up will occur; neither photolysis reactions nor spontaneous reac- 
tions can add net reducing power the atmosphere-ocean system. Our model cannot predict low steady state 
oxygen abundances in predominantly anoxic atmospheres, however; atmospheric O; is truncated at a lower 
limit of 0.1 Pa for numerical efficiency. Importantly, our modeling approach is agnostic on the plausibility 
of the various photochemical scenarios that have been proposed for abiotic O, atmosphere, such as O5- and 
CO- rich atmospheres maintained by continuous CO; photodissociation (Gao et al., 2015; Hu et al., 2020). 
Studies of these scenarios enforce global redox balance at the boundaries of the atmosphere-ocean system 
and determine whether appreciable atmospheric oxygen exists in the resultant photochemical steady state 
(e.g., Harman et al., 2015). Here, we are instead using a time-dependent model to investigate whether slight 
imbalances in atmosphere-ocean boundary fluxes can result in atmospheric oxygen accumulation on long 
timescales (cf., Luger & Barnes, 2015; Schaefer et al., 2016; Wordsworth et al., 2018). 


There are many uncertain parameterizations and parameter values in our model, and so all results are 
presented as Monte Carlo ensembles that randomly sample a wide range of uncertain parameter values. 
Parameter ranges and their justifications are described in full in the supporting information (Table S1). 
We sampled a range of temperature-dependent mantle viscosities, efficiencies of XUV-driven escape, un- 
certain early sun XUV fluxes, carbon cycle feedbacks, deep hydrological cycle dependencies, and albedo 
parameterizations. Unknown parameters that are particularly important for oxygen false positives include 
the dry crustal oxidation efficiency, faz o4, Which is the fraction of Fe** in newly produced crust that is 
oxidized to Fe** in the presence of an oxidizing atmosphere via non-aqueous reactions. This parameter is 
sampled uniformly in log space from 10 * to 10% (Gillmann et al., 2009). Another important parameter is 
the XUV-driven escape efficiency, &owxuv, which is the fraction of stellar XUV energy that drives H-escape. 
This is sampled uniformly from 0.01 to 0.3, and the portion of energy that goes into escape once the XUV 
flux exceeds what is required for O-drag is an additional free parameter. 


2.1. Model Validation 


To validate the model, we first show that it can successfully reproduce the atmospheric evolution of Earth 
and Venus. Venus results are described in detail in supporting information Section C, and here we summa- 
rize key results for Earth. Figure 2 shows Monte Carlo model outputs over a range of Earth-like volatile in- 
ventories, specifically an initial water content of 1-10 Earth oceans, an initial CO, content of 20-2,000 bar. 
Moreover, only initial CO; inventories less than the initial water inventory by mass are permitted, and 
an initial (post core-formation) mantle redox state around the Quartz-Fayalite-Magnetite (QFM) buffer is 
assumed. There is evidence for more reducing Hadean continental crust (Yang et al., 2014), and other ter- 
restrial planets such as Mars likely have reducing mantles (Wadhwa, 2001). While a rapidly oxidized mantle 
(e.g., Zahnle et al., 2010) is assumed in all nominal calculations, the sensitivity of our results to initial man- 
tle redox is explored in supporting information Section G. 
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Figure 2. Earth's coupled redox-thermal-climate evolution (without life). The model is applied to the Earth from magma ocean to present with initial water 
inventories ranging from 1 to 10 Earth oceans, and initial CO; inventories ranging from 20 to 2,000 bar. Additionally, we only plot model runs where the initial 
water inventory exceeds the initial CO; inventory. The lines are median values and shaded regions denote 95% confidence intervals across 3,000 model runs. 

In the absence of life, Earth's atmosphere after 4.5 Ga is always anoxic (c) because outgassing and crustal hydration sinks overwhelm oxygen production via 
photolysis and diffusion-limited hydrogen escape (i). (a, b) The magma ocean persists for a few million years, consistent with previous studies. (e) The magma 
ocean ends when the planet's interior cools such that heatflow from the interior drops below the runaway greenhouse limit. (d) When this occurs, liquid water 
oceans condense onto the surface, (f) a temperate carbon cycle commences. (c) There is sometimes a brief spike in atmospheric oxygen following magma ocean 
solidification due to the persistence of a steam atmosphere and hydrogen escape, (i) but this oxygen is rapidly drawn down by geological sinks. (g) Volatile 
cycling is controlled by the rate at which fresh crust is produced. Mantle redox evolution is plotted (h) alongside proxy estimates (O'Neill et al., 2018; Trail 


et al., 2011). 
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Figure 2a shows the time-evolution of mantle potential temperature and surface temperature, Figure 2b 
shows the solidification of the magma ocean from core to surface, which takes several Myr, Figure 2c shows 
the evolution of atmospheric volatile inventories, Figure 2d shows the globally averaged depth of liquid 
water oceans. The inflection in temperature evolution and solidification radius around 10* years reflects 
the transition from a low viscosity, rapidly convecting magma ocean, to more solid-like magma mush con- 
vection (Lebrun et al., 2013). Note the transfer of water from steam atmosphere (Figure 2c) to liquid water 
ocean (Figure 2d) following magma ocean solidification at around 10° years. Figure 2e shows the planetary 
energy budget, Figure 2f shows carbon outgassing and weathering fluxes, Figure 2g shows total crustal 
production, Figure 2h shows the evolution of (solid) mantle oxygen fugacity relative to the QFM buffer, and 
Figure 2i shows oxygen fluxes into/out of the atmosphere, excluding loss of oxygen to the magma ocean, 
which is modeled as instantaneous melt-atmosphere equilibration rather than a continuous sink flux; this 


temperature-dependent equilibrium partitioning controls atmospheric pO, for the first few million years 
(Figure 2c). 


Our modeled early Earth atmosphere-thermal-climate evolution is broadly consistent with semiquantita- 
tive reconstructions of Hadean atmospheric evolution (Zahnle et al., 2007, 2010). We also find that in vir- 
tually every model run, after 4.5 Gyr of atmospheric evolution, the atmosphere is anoxic (Figure 2c). This 
is unsurprising. Hydrogen escape during the initial Myr magma ocean does not add free oxygen to the 
atmosphere but instead oxidizes the interior, as has been described previously (Hamano et al., 2013). In 
some cases, small amounts of abiotic oxygen are produced in the post magma-ocean steam atmosphere, but 
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Figure 3. Oxygen false positives from high initial CO;:H;O inventories (Scenario 1). The model is applied to the Earth from magma ocean to present with 
randomly sampled initial water inventories ranging from ~0.1 to 10 Earth oceans, and initial CO; inventories ranging from ~20 to 2,000 bar (implying 
CO;:H50 ranging from 0.01 to 100 by mass). Only model outputs with modern day atmospheric oxygen exceeding 10" kg (>~0.02 bar) are plotted. Subplots 


are the same as in Figure 2, and shaded 


regions denote 95% confidence intervals. (a) High atmospheric CO; ensures the surface temperature always exceeds 


the critical point of water after the pre-main sequence, and (d) thus permanent liquid water oceans do not condense. The lack of surface water, low volatile 


content of the mantle, and high surface 


pressure increasing volatile solubility in partial melts all limits oxygen sinks. (a, i) The largest atmospheric sink is dry 


crustal oxidation, which diminishes with time as the interior cools. (c) Atmospheric oxygen produced via H escape may start to accumulate after several Gyr of 


evolution. 
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this atmospheric oxygen is rapidly overwhelmed by outgassing and other crustal sinks. Subsequent oxygen 
production via diffusion-limited escape is small, and so there are no further opportunities for abiotic accu- 
mulation so long as the planet remains geologically active. 


For Venus, the model can recover current atmospheric conditions assuming the initial water inventory is 
small, and that crustal sinks of oxygen are efficient (see supporting information Section C). Venusian histo- 
ries in which the surface was never habitable and in which the surface was habitable for several billion years 
can both be reconciled with the current atmosphere, which is broadly consistent with previous modeling of 
Venus’ atmospheric evolution (Chassefiére et al., 2012; Kasting & Pollack, 1983; Way et al., 2016). 


3. Results 


If Earth's initial volatile inventories are varied, then oxygen-rich atmospheres may be possible. Here, we 
outline three scenarios whereby an abiotic Earth could have accumulated an oxygen-rich atmosphere after 
4.5 Gyr. None of these scenarios guarantee an oxygen-rich atmosphere; instead, oxygenated atmospheres 
are a possible outcome that is dependent on the efficiency of oxygen crustal sinks and atmospheric escape. 


3.1. Scenario 1: High CO,:H,O Initial Inventory Leading to Perpetual Runaway Greenhouse 


Figure 3 shows selected model outputs for planets with initial CO;:H;O volatile inventories greater than 
one by mass and atmospheric O, > 10” kg at present (Po; > ~0.02 bar). For these planets, the greenhouse 
warming from a dense CO; atmosphere ensures that the surface temperature is above the critical point of 
water; liquid water never condenses on the surface at 1 AU, except briefly, and in small amounts, during 
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Figure 4. Conditions required for Scenario 1 oxygen false positive. Each dot denotes a single model run, and model 
runs are shown for uniformly sampled initial volatile abundances: 10-10? kg CO; and 107-10? kg HO. Atmospheric 
oxygen at 4.5 Gyr is plotted as a function of (a) dry crustal oxidation efficiency, and (b) the initial CO;:H;O inventory 
by mass. Note that Scenario 1 oxygen accumulation (high CO, perpetual runaway greenhouse atmospheres) requires 
both an initial CO2:H2O ratio >1 (green box) and for dry crustal oxidation to be relatively inefficient, with «0.195 of 
Fe** in newly produced crust oxidized. The observed range in carbonaceous chondrite CO;:H5O ratios (purple interval) 
is shown in (b) as a rough proxy for Earth's initial volatile inventory. The outliers with high oxygen (red box) are 
Scenario 3 (desertworld) false positives, which are examined in Section 3.3. Anoxic atmospheres truncate at 10 

bar for numerical efficiency (see supporting information); these model runs represent outcomes with essentially no 
atmospheric oxygen. 


early solar evolution (Figure 3d). These high CO;:H;O perpetual runaway atmospheres have been described 
previously (Marcq et al., 2017; Salvador et al., 2017). The lack of liquid surface water precludes CO;-draw- 
down via silicate weathering (Figure 3f). Reactions between supercritical water and silicates will be severely 
kinetically limited by sluggish solid state diffusion, and are therefore assumed to be negligible (Zolotov 
et al., 1997). Consequently, a dense CO; atmosphere and supercritical surface temperature persist indefi- 
nitely (Figure 3a), despite the planet residing in the habitable zone. Moreover, there is sufficient steam in 
the atmosphere to ensure diffusion-limited hydrogen escape provides an appreciable source flux of oxygen 
(Figure 3i). 


Oxygen accumulation also requires limited oxygen sinks, and this may occur on high CO;:H5O worlds for 
two reasons. First, since permanent oceans do not condense, it is difficult to sequester outgassed volatiles 
left over from the magma ocean in the interior; hydration reactions and carbonatization reactions do not oc- 
cur without liquid surface water. Limited mantle regassing after magma ocean outgassing implies low man- 
tle volatile content, which inhibits the capacity of outgassed reductants to draw down oxygen. Second, the 
high pressure from the dense CO, atmosphere, while not enough to significantly increase the silicate soli- 
dus, does increases the solubility of volatiles in partial melts. In combination with low mantle volatile con- 
tent, this ensures limited outgassing of reducing species (Figures 3f and 3i, cf., Gaillard & Scaillet [2014]). 


However, even with limited outgassing, new crust is still being produced that may be directly oxidized by 
gaseous O; (Figure 3g). Figure 4a shows atmospheric oxygen abundances after 4.5 Gyr as a function of dry 
crustal oxidation efficiency, fary—oxia, for a large number of model runs sampling 109-10? kg initial CO; and 
H30 (or ~20-2,000 bar and 0.1-10 Earth oceans). The efficiency parameter is necessarily quite low (<0.1%) 
in all the model runs that produce significant oxygen. The plausibility of such inefficient crustal oxidation is 
explored in the discussion. Figure 4b shows 4.5 Gyr atmospheric oxygen as a function of the initial CO;:H;O 
ratio, and confirms that high oxygen accumulation only occurs when CO;:H,O exceeds unity. 
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Figure 5. Oxygen false positives on waterworlds (Scenario 2). The model is applied to the Earth from magma ocean to present with randomly sampled initial 
water inventories ranging from 10 to 230 Earth oceans, and initial CO; inventories ranging from —20 to 6,000 bar. Only model outputs with modern day 
atmospheric oxygen exceeding 10" kg (>~0.02 bar) are plotted. Subplots are the same as in Figure 2, and shaded regions denote 95% confidence intervals. (d) 
The large surface volatile inventory increases the (g) mantle solidus such that melt production and tectonics shut off shortly after formation. (i) This shuts down 
all oxygen sinks and (c) allows for the gradual accumulation of oxygen via diffusion-limited H escape over several Gyr. 


3.2. Scenario 2: Waterworlds 


Figure 5 shows selected model outputs for planets with H5O volatile inventories between 10 and 230 Earth 
oceans. For these planets, a liquid water ocean condenses out of the atmosphere after a few million years 
and any oxygen left over from the post-formation steam-dominated atmosphere is typically removed by 
oxygen sinks (Figure 51), just as in the nominal Earth model (Figure 2). However, the pressure overburden 
from the large surface water inventory dramatically increases the solidus of the silicate mantle. This effect, 
which has been described previously (Kite & Ford, 2018; Noack et al., 2016), causes fresh crustal production 
to cease completely after a few billion years when the mantle potential temperature drops below the solidus 
(Figure 5g). The cessation of crustal production suppresses all geological oxygen sinks; hydration reactions 
stop as there is no melt production or new crust to oxidize (Figure 5i). The source flux of oxygen is low in 
this scenario. An effective cold trap ensures oxygen production rates of ~0.01 Tmol/yr via diffusion-limited 
escape. However, since oxygen sinks are negligible, this small source flux is sufficient to accumulate mod- 
ern Earth-like oxygen abundances over several Gyr (Figure 5c). There are also a small number of model 
runs where oxygen persists from the magma ocean, since equilibrium oxygen fugacities are high at the 
elevated solidus temperature under high pressures (Figure 5c). 


Figure 6 shows the oxygen abundance after 4.5 Gyr for individual model runs as a function of the initial wa- 
ter inventory. Waterworld oxygen false positives only begin to become likely when the initial water invento- 
ry exceeds around 50 Earth oceans or 10? kg (for Earth-sized planets). There is significant scatter in results 
due to uncertainty in the temperature-dependent mantle viscosity, which controls the duration of tectonics. 


KRISSANSEN-TOTTON ET AL. 10 of 20 


Aa l 
AJU 
ADVANCING EARTH 
AND SPACE SCIENCE 


1073 


Final O2 (bar) 


si 
i 


1077 


10? 10? 
Initial H20 (Earth oceans) 


Figure 6. Prevalence of Scenario 2 waterworld false positives. Each dot 
denotes a single model run, and model runs are shown for uniformly 
sampled initial volatile abundances: 10-107? kg CO, and 107-10% kg 
H,O. Atmospheric oxygen at 4.5 Gyr is plotted as a function of the 
initial planetary water inventory. Waterworld oxygen false positives are 
unlikely unless the initial water inventory exceeds 10? kg (~50 Earth 
oceans). There is a high probability of abiotic oxygen accumulation 

for water inventories exceeding 100 Earth oceans. Clustering around 
0.1 bar occurs because this is the amount of oxygen accumulated after 
4.5 Gyr for temperate surface conditions and Earth-like stratospheric 
water abundances. Warmer surface oceans (~400-600 K) result in more 
stratospheric water vapor and thus enable higher oxygen levels. 
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3.3. Scenario 3: Desertworlds 


The final scenario whereby abiotic oxygen could accumulate on habita- 
ble zone planets around Sun-like stars occurs for planets with extreme- 
ly small initial volatile inventories (initial water inventory < ~0.3 Earth 
oceans). Figure 7 shows selected model outputs representing this de- 
sertworld false positive. The required sequence of events are as follows: 
the low volatile inventory ensures that the magma ocean freezes quickly 
(typically ~10° years, Figure 7b), even though the planet is still in a runa- 
way greenhouse state due to high heatflow from the interior (Figure 7e). 
A steam-dominated atmosphere can therefore persist for a few million 
years, and oxygen may accumulate during this time because there is no 
surface magma ocean to dissolve the oxygen. Dry crustal oxidation will 
remove some oxygen during this steam atmosphere phase, but oxidation 
will be limited by the rate at which oxygen can diffuse into extrusive lava 
flows (Figures 7i and S1d). When a shallow ocean does eventually con- 
dense out as heatflow from the interior drops below the runaway green- 
house limit (Figures 7d and 7e), oxygen may persist for billions of years 
if oxygen sinks are small (Figure 7c). Outgassing sinks are limited by the 
low volatile inventory of the planet, but inefficient dry crustal oxidation is 
also required post-magma ocean for the oxygen to persist for 4.5 Gyr; the 
efficiency parameter, fa, 4, must be «0.196 (Figure S1). This habitable 
scenario is qualitatively different to the uninhabitable Scenario 1 false 
positives: in the former, atmospheric oxygen accumulates early and grad- 
ually declines due to crustal sinks, whereas in the latter, oxygen accumu- 
lation takes several Gyrs. Desertworld false positives also require efficient 
XUV-driven hydrodynamic escape (71090) during the steam atmosphere 
phase to produce large atmospheric oxygen abundances after the magma 
ocean has solidified (Figure S1). 


The modeling approach adopted in this paper has several important caveats and limitations. First, we con- 
sider the reasons why abiotic oxygen accumulation could be underestimated in our model. 


4.1. Assumptions That May Underestimate Abiotic Oxygen Accumulation 


As noted above, the model does not track nitrogen fluxes and instead assumes 1 bar N, partial pressure 
throughout. This limits oxygen accumulation by providing a non-condensable background gas to throttle 
hydrogen escape at the cold trap (Kleinbóhl et al., 2018; Wordsworth & Pierrehumbert, 2014). Nitrogen at- 
mospheric evolution for terrestrial planets is highly uncertain, and the evolution of Earth's atmospheric N3 
inventory is poorly constrained (Johnson & Goldblatt, 2018; Stüeken et al., 2016). However, for terrestrial 
planets that form with low nitrogen inventories, or with most of their nitrogen sequestered in the interior 
(Wordsworth, 2016), then oxygen accumulation on temperate planets could be a more common outcome 


than our modeling suggests. 


Our nominal outgassing model may overestimate fluxes of reduced gases per unit mass partial melt. This is 
because we do not account for graphite saturation and redox-dependent partitioning of carbon-bearing spe- 
cies between crystalline and melt phases; we instead assume a constant partition coefficient for relating solid 
mantle CO; content and total melt plus gas phase concentrations (e.g., Lebrun et al., 2013). Models of Mar- 
tian outgassing (Grott et al., 2011), and more generalized terrestrial outgassing models (Ortenzi et al., 2020) 
both show that reducing mantles tend to outgas fewer volatiles by mass than more oxidized mantles for 
the same amount of crustal production. Our model does, however, account for the greater reducing power 
of volcanic outgassing on planets with lower mantle oxygen fugacities (Gaillard & Scaillet, 2014; Wogan 
et al., 2020). Consequently, our conservative approach maximizes fluxes of outgassed reductants, and may 
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Figure 7. Oxygen false positives on desertworlds (Scenario 3). The model is applied to the Earth from magma ocean to present with randomly sampled 

initial water inventories ranging from ~0.05 to 0.35 Earth oceans, and initial CO; inventories ranging from ~6 to 60 bar. Only model outputs with modern day 
atmospheric oxygen exceeding 10" kg (>~0.02 bar) are plotted. Subplots are the same as in Figure 2, and shaded regions denote 95% confidence intervals. (b) 
The low initial volatile inventory ensures the magma ocean solidifies before the runaway greenhouse is over, (c) allowing for significant oxygen accumulation in 
the Myr steam atmosphere that (e) persists until the heatflow from the interior drops below the runaway greenhouse limit. If subsequent oxygen sinks are low, 
then the few bar oxygen that accumulate early on may persist for billions of years. 
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underestimate oxygen accumulation. Sensitivity tests which consider more reducing mantles and graphite 
saturated melts are described below. 


Finally, our model may underestimate the duration of steam atmospheres following magma ocean solid- 
ification, and therefore underestimate oxygen accumulation prior to ocean condensation. Based on the 
time required to precipitate an Earth ocean and the apparent absence of stable climate states at the runa- 
way greenhouse limit (Figure S4), it is typically assumed that the time required to transition from steam 
atmosphere to surface water ocean is ~10° years (Abe, 1993; Zahnle et al., 2007). In our model, the at- 
mosphere-ocean system is assumed to be in radiative equilibrium with a negligible heat capacity (Lebrun 
et al., 2013). However, for waterworlds with hundreds of Earth oceans, the time required for a steam at- 
mosphere to condense is long enough for abiotic oxygen accumulation to be significant. Moreover, stable 
climate states may exist in between a molten surface and a temperate surface ocean, especially if clouds— 
which we ignore—are included in radiative transfer calculations (Marcq et al., 2017 their Figure 6). Finally, 
in our model, oxygen is partitioned between the magma ocean and the atmosphere assuming chemical 
equilibrium. This is a reasonable assumption for high temperature/low-viscosity magma oceans with very 
short mixing times, but as the surface temperature approaches the solidus, the “magma ocean mush” will 
behave more like a solid (Lebrun et al., 2013; Salvador et al., 2017), and oxygen produced during the final 
stages of the magma ocean may not be efficiently sequestered in the melt. 


4.2. Assumptions That May Overestimate Abiotic Oxygen Accumulation 


Next, we consider model assumptions that might cause abiotic oxygen accumulation to be overestimated. 
The omission of a nitrogen cycle means the model ignores the removal of atmospheric oxygen in the ocean 
as dissolved nitrate. Nitrate formation is thermodynamically favorable at temperate surface conditions 
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(Krissansen-Totton et al., 2016), and the reaction may occur via lightning. Indeed NO-formation via light- 
ning has previously been established as an important mechanism for removing photochemically produced 
atmospheric oxygen and inhibiting oxygen false positives (Harman et al., 2018). However, nitrogen fixation 
via lighting is unlikely to prevent abiotic oxygen accumulation on waterworlds. In the absence of any other 
sources or sinks nitrate-formation via lightning could draw down the modern Earth's atmospheric oxygen 
reservoir in 20-200 Myr (Krissansen-Totton et al., 2016) and could therefore remove 1 bar of nitrogen every 
0.17-1.7 Gyr. Once most atmospheric nitrogen is converted to nitrate in the ocean, oxygen may accumulate 
rapidly due to the lack of a non-condensible cold trap (Wordsworth & Pierrehumbert, 2014), and nitrogen 
will not be significantly replenished by outgassing on waterworlds due to the pressure overburden effect 
precluding new crustal production. Waterworlds with large initial nitrogen atmospheric inventories (10 s 
of bar) could avoid abiotic oxygen accumulation if oxygen drawdown via lightning exceeds oxygen produc- 
tion via diffusion-limited hydrogen escape. However, it is debated whether nitrate is the kinetically stable 
form of nitrogen on habitable worlds (Hu & Diaz, 2019; Ranjan et al., 2019; Wong et al., 2017). Efficient 
conversion of nitrate back to molecular nitrogen via abiotic chemodentrification might return oxygen to 
the atmosphere and prevent its accumulation in the ocean, regardless of the initial N; volatile inventory. In 
summary, it is unlikely that nitrate sinks will always preclude abiotic oxygen accumulation on waterworlds, 
but better constraints on aqueous nitrogen chemistry would improve this assessment. 


Another important limitation of our model is that it ignores infrared cooling of the upper atmosphere and 
the throttling of hydrogen escape that results from a cool stratosphere. In our nominal model, an isothermal 
210 K stratosphere is assumed. However, CO;-rich atmospheres may efficiently radiate in the IR, cooling the 
stratosphere and enhancing the water cold trap (Wordsworth & Pierrehumbert, 2013). To test the sensitivity 
of our results to stratospheric temperature, we repeated our calculations and introduced an additional strat- 
ospheric temperature variable, which was randomly sampled from 150 to 250 K. The full results of these 
calculations are shown in Figure S13. To summarize, neglecting the stratospheric radiation budget does 
not affect the viability of Scenarios 2 (waterworlds) or 3 (desertworlds) because, in the former, the oxygen 
source is diffusion-limited escape through a N;-dominated atmosphere at modern Earth-like rates, whereas 
for the later, oxygen accumulation occurs predominantly during the early, steam-dominated atmosphere, 
and so stratospheric temperature does not have a strong influence on escape rates (Figures S13b and S13c). 
However, oxygen accumulation via Scenario 1 (high CO;:H50 perpetual runaway greenhouse) is unlikely if 
the stratosphere is cooler than 200 K (Figure S132). Photochemically produced ozone or hazes may offset 
the cooling effect of CO;, and so a full assessment of Scenario 1 requires more detailed radiative-photo- 
chemical modeling. 


One additional caveat is that our model assumes an anhydrous solidus. This simplification is probably rea- 
sonable for post-magma ocean mantle conditions (Kite & Ford, 2018), but it is possible to imagine a scenario 
whereby waterworld mantles become increasingly hydrated via subduction after the magma ocean phase, 
and that this hydration offsets the pressure overburden effect to maintain geologic activity, and therefore 
oxygen sinks, for much longer than our nominal model suggests. To test this possibility, we conducted a sen- 
sitivity test where we accounted for mantle hydration decreasing the solidus (Katz et al., 2003). The results 
of this sensitivity test are discussed in detail in supporting information Section D, but in summary, mantle 
hydration does not have a large effect on oxygen accumulation on waterworlds; it merely shifts the ocean 
mass threshold for oxygen accumulation. 


Sensitivity tests were also conducted to assess whether the delivery of reducing material such as metallic 
iron and FeO via impacts (Zahnle et al., 2020), a more reducing initial mantle, or larger planet-star sep- 
arations could inhibit abiotic oxygen accumulation. These sensitivity test results are described in full in 
supporting information Sections E, G, and H, respectively. In summary, we find that high impactor fluxes 
could preclude desertworld oxygen accumulation assuming all impactor material is completely oxidized 
(supporting information Section E). Impactors may thus prevent Scenario 3 (desertworld) false positives in 
some cases, although this is not guaranteed because it is possible to imagine planetary formation pathways 
with smaller impactor fluxes and/or where the majority of impactor material is buried or lost to space as 
suggested by some impact simulations (Marchi et al., 2018). Scenarios 1 (high CO2:H2O perpetual runa- 
way greenhouse) and 2 (waterworlds) are viable under a more reducing (iron-wüstite buffer) initial man- 
tle (supporting information Section G). This counterintuitive result occurs because, even though degassed 
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volatiles are likely to be more reducing, total volatile concentrations in the melt phase are typically lower 
due to graphite saturation (Grott et al., 2011; Ortenzi et al., 2020). Moreover, crustal sinks are precluded 
by high overburden pressure, regardless of the redox state of the crustal material. Although Scenario 3 
(desertworlds) is seemingly not excluded by a more reducing mantle, evaluating this would require more 
complete radiative transfer and photochemical modeling of CO-H, dominated atmospheres. Finally, when 
nominal calculations are repeated at 1.3 AU, both Scenario 2 and Scenario 3 oxygen false positives still occur 
frequently (supporting information Section H). Scenario 1 false positives do not occur at large planet-star 
separations because a high CO,:H,O atmosphere cannot maintain a perpetual runaway greenhouse state 
after magma ocean solidification. 


For Scenario 1 and 3 to be viable, the dry crustal oxidation parameter, fa; 444, must be relatively small 
(<0.1%). This contrasts with Venus where fa oia probably needs to exceed >0.1% to remove virtually all 
O, from the atmosphere (see Venus validation in supporting information). This parameter is challenging to 
definitively constrain because it represents a range of physical processes including the diffusion of oxygen 
into extrusive lava flows (Gillmann et al., 2009), direct oxidation of small grain erosion products (Arvidson 
et al., 1992), and various other gas-solid redox reactions (Zolotov, 2019). Even if the oxidation of fresh crust 
is typically efficient, low dry crustal oxidation efficiencies cannot be ruled out because tectonic regimes 
where most magmatic activity is intrusive and isolated from the atmosphere are possible. In any case, the 
uncertainties in crustal oxidation processes highlights need for future missions to Venus to better constrain 
its redox evolution. Mars' crust is more oxidized than its upper mantle, but this oxidation cannot necessarily 
be used to constrain dry crustal oxidation efficiency since it might be attributable to early hydrous alteration 
(Herd et al., 2002; Wadhwa, 2001). The presence of gray, reduced sediments mere centimeters below more 
oxidized Martian regolith, as revealed by Curiosity, argues against efficient post depositional gas-solid oxi- 
dation under oxic conditions (Ming et al., 2014). 


Finally, we note that in some cases abiotic oxygen accumulation is contingent on highly uncertain atmos- 
pheric escape physics. This uncertainty does not affect the viability of the waterworld (Scenario 2) false 
positives because the required H escape flux is comparable to the modern Earth's diffusion-limited escape 
flux (Catling & Kasting, 2017, p. 148). Oxygen accumulation only occurs in this case because of the pres- 
sure overburden suppression of oxygen sinks. Scenario 1 (high CO;:H50 perpetual runaway greenhouse) 
is similarly unaffected. However, for Scenario 3 (desertworlds), oxygen accumulation only occurs because 
of efficient XUV-driven escape of hydrogen, £jyxuy > 0.1. If H-escape is photochemically limited (e.g., 
Wordsworth et al., 2018) then oxygen accumulation may be limited by the photochemical dissociation of 
water by UV photons, and the rate and which H5O recombination reactions occur. However, the oxygen 
source fluxes required in our desertworld scenario (~500 Tmol O;/yr, Figure 7i) are comparable to the 
water loss rates inferred for a steam-only early Earth atmosphere calculation using a photochemical model 
(Wordsworth et al., 2018, their Figure 9). Future work ought to couple the geochemical evolution model 
here to a photochemical model that includes C-bearing species to better assess the potential for oxygen 
accumulation on desertworlds. Finally, it is possible that non-thermal O loss (Airapetian et al., 2017) or pho- 
tochemically modulated stoichiometric escape of H and O (McElroy, 1972) could lessen oxygen accumula- 
tion. Upcoming JWST observations of highly irradiated terrestrial planets may constrain escape processes 
and improve predictions of oxygen accumulation on more temperate planets. 


Sulfur outgassing and burial may have played an important role in the oxygenation of Earth's atmosphere 
(Gaillard et al., 2011; Olson et al., 2019). Sulfur species are ignored in our nominal model because their 
bulk abundances are probably too small to qualitatively change our oxygenation scenarios (Wordsworth 
et al., 2018). Supporting information Section I explores the consequences of adding reduced sulfur species 
to our outgassing model, and confirms that, for Earth-like sulfur mantle abundances, total oxygen sinks are 
comparable to when sulfur is neglected. With that said, mantle sulfur abundances are contingent on forma- 
tion processes (e.g., Grewal et al., 2019) and could be highly variable. Incorporating a complete sulfur cycle 
into a redox evolution model to investigate the sensitivity of oxygenation to initial sulfur content is an op- 
portunity for future research. Note, however, that mantle sulfur abundances are irrelevant for waterworlds 
(Scenario 2) where all crustal production is suppressed. 


In summary, there are several unknowns that preclude definitive predictions of how frequently the three 
scenarios outlined in this study might occur, but none can be ruled-out with current knowledge. 
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4.3. Implications for Future Observations 


How might future observations discriminate between the three abiotic oxygen scenarios described above 
and oxygen produced by a biosphere? In principle, high CO;:H5O atmospheres should be possible to diag- 
nose via direct imaging spectral observations because they are not habitable. A clear atmosphere is likely 
since the coexistence of atmospheric H5O, O2, O3, and abundant OH radicals may preclude the accumula- 
tion of photochemical hazes. Strong CO; absorption features ought to be visible, as should pressure-sensi- 
tive CIA features from the high pressures; more detailed photochemical and spectroscopic simulations will 
be required to determine the best false positive discriminants for these worlds. 


Because waterworlds and desertworlds are habitable, they may be more challenging to discriminate from 
inhabited terrestrial planets. Crucially, waterworld false positives would be ruled out by a detection of sub- 
aerial land because, for Earth-like gravity, the presence of emerged continents limits the maximum ocean 
depth to around 10 km (Cowan & Abbot, 2014), or equivalently, a few Earth oceans by mass. This limit 
arises because silicates cannot support their own weight with greater topography. Consequently, the detec- 
tion of an ocean-continent dichotomy using time-resolved photometric mapping (Cowan et al., 2009; Farr 
et al., 2018; Fujii et al., 2010; Kawahara & Fujii, 2010; Lustig-Yaeger et al., 2018) could rule out a waterworld 
false positive, assuming alternative explanations for dichotomies in surface maps could be excluded. This 
highlights the need for large aperture direct imaging mission to ensure sufficient time-resolution to map the 
surface over a planet's rotation. Alternatively, independent mass and radius constraints from radial velocity 
observations and thermal infrared direct imaging (Quanz et al., 2019), respectively, could also help rule out 
large (few wt.%) water inventories based on bulk density. 


Desertworlds are likely the most challenging scenario to disambiguate from biological oxygen. Time-re- 
solved photometric surface maps and/or the lack of ocean glint could help evaluate the surface water inven- 
tory and might be suggestive of a small water inventory (Lustig-Yaeger et al., 2018; Robinson et al., 2010). 
There are potentially other diagnostic spectral signatures of desertworlds such as spatial variation in at- 
mospheric water vapor and photochemistry that could be tested using general circulation models and pho- 
tochemical models. The presence of long-lived sulfuric acid hazes (Loftus et al., 2019) has been proposed 
as putting an upper bound on surface water abundances, but the desertworlds considered here likely have 
larger surface water inventories than this threshold. 


Broadly speaking, the scenarios outlined in this study emphasize that no single observation, including oxy- 
gen detection on habitable zone planets around sun-like stars, will be uniquely diagnostic of life. It will be 
necessary to design future telescopes that are capable of both constraining the full planetary/stellar context 
and identifying multiple lines of evidence for life (Catling et al., 2018; Walker et al., 2018). For example, 
oxygen detection on an ostensibly habitable terrestrial planet would be persuasive if accompanied by sur- 
face biosignature detections (Schwieterman et al., 2018), temporal biosignatures (Olson et al., 2018), or 
co-existing reducing gases in atmospheric disequilibrium (Krissansen-Totton et al., 2016). The coexistence 
of oxygen and methane remains an excellent biosignature and would not be expected for any of the oxygen 
false-positive scenarios described above. Indeed, it is difficult to produce large methane abundances in hab- 
itable planet atmospheres without life, even in anoxic atmospheres (Krissansen-Totton, Olson, et al., 2018; 
Wogan et al., 2020). It should also be noted that the scenarios in this study were illustrated for habitable zone 
planets around sun-like stars, but they may also be applicable to habitable zone planets around M-dwarfs. 


5. Conclusions 


The redox evolution of habitable zone terrestrial planets is strongly dependent on initial volatile inven- 
tories and the efficiency of crustal sinks. Uninhabited, Earth-sized planets within the habitable zone of 
G-type stars are very unlikely to accumulate abiotic oxygen if their initial volatile inventories are Earth- 
like. However, if initial volatile inventories differ dramatically from that of the Earth, then non-biological 
oxygen accumulation is possible, even when atmospheric noncondensable inventories are large. This may 
occur when either (i) the initial CO;:H50 ratio exceeds one, which suppresses oxygen sinks due to the low 
mantle volatile content and because surface conditions are too hot for aqueous reactions, or (ii) the initial 
H20 inventory is very large, thereby halting crustal production after a few billion years and shutting off all 
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oxygen sinks, or (iii) the planet is very volatile-poor, in which case oxygen may accumulate during the steam 
atmosphere that persists after magma ocean solidification. Inefficient dry crustal oxidation is required for 
scenarios (i) and (iii) to yield large oxygen abundances, and scenario (i) is sensitive to stratospheric temper- 
ature. Fortunately, observational discriminants exist for all three of these scenarios; scenario (i) planets are 
uninhabitable, whereas the ability to constrain surface water inventories using time-resolved photometry 
would be useful for ruling out scenarios (ii) and (iii). More generally, the possible existence of these oxy- 
gen false positive scenarios highlights the need for a systems approach to biosignature assessment where 
biogenicity is judged not by the presence or absence of a single biosignature gas, but by multiple lines of 
evidence from both spectrally resolved and temporally resolved observations. 
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Fig. S1. Conditions required for Scenario 3 oxygen false positive. Each dot denotes a 
single model run, and model runs are shown for all model runs in Fig. 7. Atmospheric 
oxygen at 4.5 Gyrs is plotted as a function of (a) dry oxidation efficiency, and (b) XUV- 
driven escape efficiency. Subplot (c) shows the initial water and carbon dioxide 
inventories that result in desertworld false positives, and (d) shows the maximum 
fractional planetary area that is continuously molten (available for oxidation) after the 
magma ocean has solidified. 
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Fig. S2: Illustrative model outputs. Selected model outputs from a single, 
representative model run taken from Fig. 2 in the main text, the nominal Earth model. 
Subplots illustrate (a) free oxygen reservoirs, (b) water reservoirs, (c) carbon dioxide 
reservoirs, (d) free oxygen fluxes, (e) atmosphere-ocean partitioning of water, (f) carbon 
outgassing and weathering fluxes, (g) iron oxidation in the solid mantle, (h) surface and 
mantle potential temperature, and (i) the solidification radius. Early fluctuations in 
surface temperature and volatile fluxes are due to non-monotonic evolution of solar 
luminosity during the pre-main sequence. Total free oxygen (a) increases and total water 
decreases (b) as H is preferentially lost to space, whereas total CO» remains constant (c). 


Text A. Model description 


A.1) Thermal evolution: 
Planetary thermal evolution is specified by energy budget and temperature-dependent 
viscosity. The time-evolution of mantle potential temperature, T, (K), is determined by 
the following equations, representing the magma ocean and solid-state convection 
phases, respectively: 
ATP nQ radioactive (r, A E rè) 
3 


l r? dr, ES A7, C, (r; -n) dT, 
fusion’ s dt 3 dt 


4ng, r; + O ors + 4ip,,H 
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ATP „Q radioactive (r? E rè) Azq, r. 2 Oy = 470,,C p (i = a) dT, (2) 
m'p core 


3 3 dt 


Here, Q 


radioactive 


(J/kg/s) is the heat production via radiogenic isotopes. Earth's radiogenic 


inventories are taken from Lebrun et al. (2013), and a range of 0.33-3.0 times this 
nominal, Earth-like inventory is sampled in Monte Carlo calculations. Additionally, p„ = 


4000 kg/m? is the assumed average density of mantle material, r, is the total radius of 
fusion B 4x 1 0? 
J/kg is the latent heat of fusion of silicates, c, = 1200 J/kg/K is the heat capacity of 


the planet, r, is the radius of the core, r, is the radius of solidification, H 


silicates. We adopt a simple exponential decay for core heatflow (J/s): 
Q „e - 20x10" x exp(-0.2x (t — 4.5 Gyrs)) (3) 


Note that we do not account for uncertainty in core heatflow since we are already 
sampling an order of magnitude range in radiogenic inventories (see above). Moreover, 
by choosing a core heatflow history at the higher end of literature estimates (Nimmo 
2007; O'Rourke & Stevenson 2016), we are effectively maximizing crustal recycling and 
subsequent oxygen sinks. Tidal heating is ignored; if Earth-moon system tidal heating 
were included then the duration of the magma ocean could be extended by a few million 
years, potentially providing more time for oxidation via H escape (Zahnle et al. 2015). 
Equation (1) governs thermal evolution except in rare cases where transition to runaway 
greenhouse causes surface temperature to increase above mantle potential temperature, 
in which case a conduction regime is adopted (see below). 


The heatflow from the convecting interior, q, (W/m?), is parameterized as follows: 


i. LJ Fag) y 


convect 


3 
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Here, æ 22x10? K' is the thermal expansion coefficient for silicates, g (m/s?) is surface 
con = 1100 is the critical Rayleigh number, k = 4.2 W/m/K is the thermal 
conductivity of silicates, 8 -1/3 and x 210? m?/s is thermal diffusivity (Lebrun et al. 
2013; Schaefer et al. 2016). Heatflow is dependent on the kinematic viscosity, v , which is 
a function of mantle potential temperature (see below). This parameterization is adopted 
for both solid state and magma ocean phases. 


gravity, Ra 


A.2) Solidus parameterization and magma ocean freezing: 


The solidus controls both the freezing of the magma ocean and the production of partial 
melt (and outgassing) during temperate geochemical evolution. The adiabatic mantle 
temperature profile, solidus, and liquidus are parameterized as follows. 


= EN NEN 
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| T(r)exp(10°(-r, +r +10°)) +7, (r)exp(10°(r, -r -107)) 
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overburden 
Tiu (r, Poverburden) = T ondes (r, Poverburden ) + 600 


Here, T, is a linear fit to the solidus for low pressure dry peridotite, and T, is a linear fit 


to the solidus for the high pressure lower mantle (Hirschmann 2000; Schaefer et al. 
2016). A smooth function between them is assumed for T. so that an analytic 


solidus 
derivative exists at all radii (see below). Following Schaefer et al. (2016), the liquidus is 
assumed to be 600 K warmer than the solidus at all pressures. We also allow for 
modulation of the solidus and liquidus by the pressure overburden of surface volatiles. 
Here, P. is the pressure from all H20 (liquid and gaseous), CO», and O» at the 


overburden 
surface. The pressure overburden is only accounted for after the magma ocean has 
solidified and after the mantle has degassed. 


The time evolution of the solidification radius is determined by a similar method to 
Schaefer et al. (2016). The rate of change in the solidification radius can be obtained by 
noting that the time derivative of solidus evolution and the time derivative of the 
adiabatic temperature profile must be equal at the solidification radius: 
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Equation (6) can be expanded using the chain rule and rearranged to obtain an 
expression for the solidus evolution in terms of the potential temperature time- 
derivative: 
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Note that the solidification radius must remain constant when the core-mantle boundary 
temperature exceeds the solidus temperature and when the surface temperature drops 
below the solidus to ensure the solidification radius and mantle potential temperature 
begin evolving together when T(r.) 2 T,,,,,(r.) and stop evolving when the magma 


ocean freezes. This procedure was checked against simply numerically solving for the 
solidification radius at every timestep: 
T olidus (r) = T(r) 


T, (r.)exp(10? (=r, tr, 105) - T, G)exp(10?(r, qii -10°)) =T Lea S (r r, ) 
exp(10? (-r, tr, +10°))+exp(10°(r, DE -10°)) B j 
(8) 


This approach is more computationally expensive but yields an identical solidification 
radius evolution to the analytic expression eq. (7). 


A.3) Mantle Viscosity Parameterization: 


Our viscosity parameterization needs to have several properties. First, it must successfully 
reproduce the modern Earth's heatflow, melt production, and plate velocity. Second, it 
needs to transition smoothly from low viscosity magma ocean, to magma mush, to solid 
state convection. Our parameterization is a variation of those assumed in other magma 
ocean-to-solid interior evolution models (Lebrun et al. 2013; Salvador et al. 2017; 
Schaefer et al. 2016). However, the parameterizations in these studies needed 
modification because they predict a low viscosity magma-ocean or mush at the modern 


Earth's potential temperature (~ 1620 K). Noting that there is a large uncertainty in the 
critical melt fraction that controls the transition from solid-like to fluid-like convection 
(Costa et al. 2009), we adopted a parameterization that ensures this transition occurs at a 
temperature that exceeds the modern Earth's potential temperature. This is illustrated in 
Fig. S3, which compares our viscosity parameterization to others in the literature. 
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Fig. S3: Assumed viscosity parameterizations compared to other parameterizations from 
the literature (solid lines). The dashed-cyan lines represent the mantle potential 
temperature and viscosity required to reproduce the modern Earth's melt production 
and plate velocity. 
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Here, V. 


coef 


= 10! to 10? Pa s is a randomly sampled parameter that accounts for 


uncertainty in solid-state viscosity. 


A.4) Stellar evolution 
We assume solar bolometric luminosity evolution, L(t) (W), for all model runs (Baraffe et 
al. 1998; Baraffe et al. 2002). For the evolution of stellar XUV luminosity, we follow the 
empirical fit developed in Tu et al. (2015). The early sun's rotation rate, Q,, is an 
unknown parameter sampled uniformly from 1.8 to 45 (relative to modern) in log space. 
From the early sun's rotation rate, the time for the early sun to fall out of saturation, t, 
(Myr), is given by: 

[929055 (10) 


For the chosen range of rotation rates, sampled saturation times range from 6 to 226 
Myrs. We assume that at saturation the sun's XUV luminosity is 10°" L(t) . To retrieve 
the modern XUV flux, we define the exponent, # —0.86/(0.351og,,(Q,) — 0.98) , and 
specify the XUV evolution as follows: 


10^? LQ), t «t, 
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The flux received at each planet's orbital distance, D (m), is calculated using 


planet—star 
Earth-sun and Venus-sun separations. We begin our model runs at ¢ = 10 Myrs in the 
stellar evolution model, but results are insensitive to the choice of zero point. 


A.5) Surface Energy Budget: 

At each time-step in the model, the surface temperature is solved numerically by finding 
the surface temperature that ensures heatflow from the interior plus absorbed shortwave 
radiation (ASR) is exactly balanced by outgoing longwave radiation (OLR). This 
assumption ignores the intrinsic heat capacity of the ocean, which is reasonable for 
Earth-like oceans, but for waterworlds this radiative equilibrium approach may 
underestimate the transition time from runaway greenhouse to surface ocean (see 
Discussion). Surface temperature is found by solving the following equation for 7,,,,: 


An (7, oT, 


surf ? 


v(T, )) + ASR(a;, s, L(t) = OLR(T,,,, T M yw no M pua. co, > [C0; ]) 
(12) 


Here, g,, is specified by equation (4). Absorbed shortwave radiation is a function of the 
planetary albedo and stellar luminosity: 

L()(1- 0,4) 
16zD A 


planet —star 


ASR(a, ui L(t)) = (1 3) 


Here, D (m) is the distance between Earth (or Venus) and the sun. Following 


planet—star 


Pluriel et al. (2019) we assumed the following parameterization for albedo: 


5 


transition ^— 


H,OY 
1000.0 + 200 x log, (ES , pH,0 » 105 Pa 


1000.0, pH,O x10? Pa (14) 
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Here, T incision (K) is the albedo transition temperature that controls the transition from a 


hot-state albedo, A,, to a cold-state albedo, A,,,,. The "hot" and “cold” states do not 
refer to non-glaciated and glaciated states, which we do not consider in our model. 
Instead, they allow for a transition from low albedo cloud free runaway greenhouse 
atmospheres at high (>~ 1000 K) temperatures, to a range of cloudy and non-cloudy 
states under more temperate conditions (Pluriel et al. 2019). This distinction is important 
for modeling Venus, where the “cold” state albedo is ~0.7, but the albedo during the 
initial runaway greenhouse phase was potentially much lower. Albedos for hot and cold 
states are sampled uniformly form 0-0.3 and 0.25-0.35, respectively for Earth and 0-0.3 
and 0.2-0.7 for Venus. The albedo of the hot state must always be equal to or less than 
that of the cold state. 


To calculate the outgoing longwave radiation (OLR), we used the publicly available code 
from Marcq et al. (2017). This code uses DISORT (Stamnes et al. 1988) with four stream 
longwave radiative transfer. The radiative transfer model only considers opacity due to 
water vapor and carbon dioxide. Rock vapor opacities are ignored since the time spent at 
rock-vaporizing temperatures is very short and unlikely to affect long term redox 
evolution. Correlated k coefficients are calculated from the high resolution molecular 
absorption spectra computed with kspectrum (Eymet et al. 2016), H2O-H;0 continuum 
absorption is taken from Clough et al. (2005), and CO2-COz continuum absorption from 
fits to Venus observations (Bézard et al. 2011). H20-CO» continuum opacity is not 
considered, and is likely negligible compared to H20-H;O and CO?-CO» continuum 
absorption (Ma & Tipping 1992). The runaway greenhouse limit calculated using the 
code of Marcq et al. (2017) closely agrees with line-by-line calculations in Goldblatt et al. 
(2013). 


The atmosphere model calculates atmospheric structure and abundance profiles using 
the expressions for dry and moist adiabats in Kasting (1988), and the thermodynamic 
properties of water are taken from steam tables (Haar et al. 1984). Given an assumed 
surface temperature and total water inventory, the code calculates the atmospheric 
water vapor profile assuming a dry convective regime (partial pressure of water vapor 
less than saturation) to moist convective regime (partial pressure of water vapor equals 
saturation) to isothermal temperature structure, where the isotherm temperature is the 
planetary skin temperature. The dry convective regime may not be present if water vapor 


is saturated at the surface. Once the water vapor profile has been calculated, the 
remainder of the surface water inventory (if any) resides in a surface water ocean. 


If a portion of the surface water resides in an ocean, then the partitioning of carbon 
between the atmosphere and ocean also determines the OLR. Thus, OLR is a function of 
dissolved carbonate concentrations. This is calculated as follows: 


[coz ]= Ox Ky Tur) (15) 


[c^ ] 


Here, rather than explicitly track the ocean alkalinity budget, we follow the approach of 
Schwieterman et al. (2019) and assume that carbonate precipitation will ensure the long- 
term carbonate saturation state of the ocean, Q, is constant. Values for Q are sampled 
randomly from 1 to 10 to allow for abiotic supersaturation. Similarly, rather than 
explicitly track cation weathering budgets, we assumed constant dissolved calcium 
abundances and sample uniformly in log space from 10^ to 3x10 mol/kg. Dissolved 
calcium concentrations in Earth oceans have varied from 10? to 3x 10'! mol/kg over 
Earth history (Halevy & Bachan 2017), but we adopt a broader range to account for 
different crustal compositions and ocean volumes. For example, Kite and Ford (2018) 
consider waterworld scenarios with essentially zero Ca^* up to 0.25 mol/kg based on 
thermodynamic models of basalt-water interaction. Explicitly tracking the ocean alkalinity 
budget is an opportunity for future research. The temperature dependent solubility 


product, K,,, is the same as in Krissansen-Totton et al. (2018). Once the total carbon 


reservoir, the ocean size, and the dissolved carbonate concentrations are known, the 
entire carbonate equilibrium system of equations can be solved to determine 
atmospheric pCO» (Krissansen-Totton et al. 2018). 


Rather than call the atmospheric radiative transfer code in real time, we precomputed a 
grid of OLR values as a function of surface temperature (250-4000 K), surface water (10 
Pa — 1 GPa), surface carbon dioxide (10 Pa to 0.1 GPa), and planetary effective 
temperature (150 — 350 K). Within the grid we linearly interpolate between grid points, 
and on the rare occasion when the model moves beyond the grid, linear extrapolation is 
adopted. A 1 bar partial pressure of N2 is assumed at every grid point. Since the 
atmospheric model calculates atmospheric structure, stratospheric mixing ratios are 
obtained, which are used to determine atmospheric escape rates (see below). 
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Fig. S4: Illustrative figure showing calculations from OLR grid. Top plots show OLR as a 
function of surface temperature for varying surface CO» pressures (log:o(Pa)). Here, was 
have assumed [CO3*] = 10? mol/kg. Middle plots show the amount of water in the 
atmosphere as a fraction of the total surface water inventory, fr. o. Bottom plots 


atmo—H 
show the stratospheric water mixing ratio for different CO» inventories. On the left-hand 


side 1 Earth ocean is assumed, whereas on the right-hand side 0.1 Earth oceans are 
assumed. 


A.6) Volatile reservoirs and planetary redox budget: 
The time evolution of planetary volatile budgets and redox states is determined by the 
following system of equations: 
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Here, M sa-i (kg) and M 
interior and the fluid magma ocean plus surface reservoirs, respectively. Water, free 
oxygen, solid FeO and FeO:; in the mantle, and carbon dioxide are separately tracked in 
the model. The variables fr represents the mass fraction of the i-th species in the 
magma ocean partial melt. These are calculated at every timestep using the equilibrium 
relations described below. Key constants include the density of mantle material, o, = 


(kg) represent the masses of the i-th species in the solid 


Fluid—i 


4000 kg/m?, molecular masses of key species, u, (kg/mol), and the assumed partition 
coefficients for CO? and H20, kyo = 0.01 and kco, = 2x10”, respectively (Lebrun et al. 


2013). See the Discussion section in the main text for possible consequences of assuming 
constant partition coefficients. Escape fluxes for atomic hydrogen, E, (mol/s’), atomic 


oxygen E, (mol/s?), and carbon dioxide, Eco, (mol/s?) are parameterized below. The 


remaining fluxes (kg/s) are the ingassing of surface water into the interior F 


ingas—H,O-gain ! 


and the loss of surface water F. 


outgassing of water from the interior, F ings Hoss 


outgas—-H5O ! 


(in general F,,, 45 4, $ F because serpentinizing reactions can remove water 


ingas- H,0— gain 
from the surface that never reaches the interior because the hydrogen produced is lost 
to space), oxidation of the interior, Fi) sig, corresponding loss of oxidants from surface 


reservoirs, F,;, ,,;, (typically F,;, ,,,, = F,;, ,,4, except for anoxic atmospheres where 
reductants are lost to space), outgassing of carbon to the surface, Fog co, and loss of 
surface carbon to weathering, F,,,,,, «.,. Fluxes of ferrous and ferric iron are equal to 

2 
E oxid —solid 
zero during solid-state convection and are described in detail in their corresponding 
sections below. Total fluid volatile masses are converted to partial pressures using 


scaled by appropriate molecular masses. All fluxes denoted F, are only non- 
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atmospheric mean molecular weight (no magma ocean) or using the melt solubility 
relationships described below (magma ocean). 


Fig. S2 shows illustrative outputs from a single model run. Subplots Fig. S2a, S2b, and 
S2c show the evolution free oxygen, water, and carbon dioxide reservoirs, respectively, as 
governed by equation (16). 


A.7) Magma-ocean evolution: 


Whilst the magma ocean exists, volatiles in fluid phases are partitioned between the 
melt, melt crystals, and the atmosphere. For water, this partitioning is described by the 
following equation (Schaefer et al. 2016): 


fio Agr, fey 1/0.74 E 
ky o fru oM orsta + (Mus ES M sut ) Frio + i g 344x103 x fluid -H50 
(17) 


Here, M puia - Ap, [3(r; -r ) (kg) is the mass of the magma ocean, M „sa - (17 V)M,, is 


the crystal mass fraction in the magma ocean, which depends on melt fraction, y (see 
below). The molecular mass of the atmosphere is given by ZZ (mol/kg). The third term on 


the left-hand side represents the mass of water in the atmosphere, where we adopt the 
solubility relationship from Papale (1997). An analogous expression can be used to 
calculate the partitioning of carbon dioxide: 
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Here, we use the solubility relationship from Pan et al. (1991). 


For oxygen, equilibrium partitioning is more complicated because both Fe** and Fe?* 
melt phases must be included. We adopt the experimental fit in Kress and Carmichael 
(1991): 
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(19) 


Here, x,,,, and x, are the mole fractions of iron-bearing species in the melt, whereas 
Xy, 22x, 9, + Xreo Is the total mole fraction of all iron-bearing species in the mantle 


(constant). The molar abundances of other species, x,, are assumed to be represent Bulk 
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Silicate Earth (White 2013). We calculate oxygen solubility using surface temperature and 
pressure conditions because the interface controls volatile partitioning. 


In equation (19) there are three unknowns Coe and Xpe o, ) and so it is necessary to 


FeO ! 
make appropriate substitutions to solve for x,,.. This can be done by expressing the 


partial pressure of oxygen in terms of the total fluid free oxygen minus the free oxygen 
dissolved in the magma ocean, and by noting that xpo = x,, - 2x, : 
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Here, the average molecular weight of silicates is taken to be: 


Hsia = XucoHmgo + Xsio, Hsio, + Xano, Hano, + Xcaoklcao + Xre,0,HFe,0, + Xreo!reo (21) 


Equation (20) is solved at every timestep in the model, and from the values for melt 
concentrations of iron-bearing species, we can retrieve the melt fractions (by mass) that 
dictate the time-evolution of volatile reservoirs: 
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Equation (19) is also used to calculate mantle oxygen fugacity for the purposes of 
outgassing calculations by substituting the solid mantle molar fractions of oxidized and 
reduced iron. 


A.8) Transition from magma ocean to solid mantle convection: 
The model switches between magma ocean and solid-state convection freely, as dictated 
by the radiation and interior heating budget. At each time step, surface temperature is 
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compared to the solidus. For as long as the surface temperature exceeds the solidus, the 
magma ocean model is adopted (equation (1)), and the solidification radius evolves with 
time according to equation (7). However, once the surface temperature drops below the 
solidus, the solidification radius is set to the planetary radius, and solid state interior 
evolution is dictated by equation (2). Volatiles are instantaneously exchanged between 
the magma ocean and the solid interior at this transition. The model assumes that when 
the surface freezes, any volatiles still dissolved in the magma mush remain in the interior 
(e.g. as basaltic glass, or gas bubbles in melt inclusions), which is reasonable given the 
short timescale for magma ocean solidification and the high viscosity of the late-stage 
magma mush. This assumption also maximizes the mantle’s capacity for subsequent 
outgassing of reduced products that remove oxygen from the atmosphere. 


During the transition from magma ocean to solid-state convection, volatile inventories 
undergo the following one-off adjustment: 
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Note that the water transferred to the mantle cannot exceed the maximum water 
content of the mantle. Transition from solid state convection to magma ocean, which is 
typically only relevant for Venus model runs, is similarly calculated.: 

M yu i "E M sui i zx Mi, M sui i [M nanie 

M nui = M mia-i E M M suia i [M nante 


liq 


(23) 


Fig. S5 shows the fraction of total CO» and H20 that reside in the solid mantle 
immediately after magma ocean solidification for outputs from Fig. 2 in the main text. 
These mantle fractions are determined by the partitioning of volatiles in the magma 
ocean as described in Supplementary Section A.7, and the instantaneous retention of 
leftover melt when surface temperature drops below the solidus, as described in this 
section. While we do not explicitly model mechanisms of volatile retention in the magma 
ocean, such as compaction within the moving freezing front, our spread of final mantle 
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volatile fractions is comparable to that of more detailed models (e.g. Hier-Majumder & 
Hirschmann 2017). 
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Fig. S5: Volatile fraction in the solid mantle immediately after magma ocean solidification 
for nominal model calculations (Fig. 2). The left subplot shows the carbon dioxide mantle 
fraction, whereas the right subplot shows the water mantle fraction. 


In a few rare cases, the transition back from solid to magma ocean causes the surface 

temperature to exceed the mantle potential temperature. When this occurs, we modify 

the energy budget as follows: 
ASR(a,,,,. L(t)) = OLR(, 


surf ? 


T, ? M rw no ? M nii co, d [CO p + qe (T, 7 Tur ) (24) 


Here, q, (W/m^?) is the conduction of heat from the surface to the interior and is 
approximated by the diffusion equation: 


Tir T T, 
lly Tang) =k E (aa) 
Pp c 


The time evolution of volatile reservoirs is also modified in this regime to account for the 
fact that the radius of solidification is moving downward towards the core (see source 
code for full details). This lasts until surface temperature is less than mantle potential 
temperature, and the model switches back to a convective mantle. 


A.9) Atmospheric Escape Parameterization: 

Atmospheric escape controls the source flux of abiotic oxygen. Escape rates are 
determined by the composition and temperature of the stratosphere, and by the stellar 
XUV flux. For low stratospheric water abundances, escape is limited by the diffusion of 
water through background gases, or by the XUV flux from the star (whatever is smaller). 
As the water content in the upper atmosphere increases, the escape regime transitions 
to XUV-limited because, for steam-dominated atmosphere, there is no cold trap limiting 
the supply of water to the upper atmosphere. Our approach does not rigorously capture 
the complexities of atmospheric escape physics (e.g. Owen 2019), but instead uses 
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plausible parameterizations that incorporate broad parameter ranges to cover a wide 
range of uncertain physical processes. Moreover, our parameterization collapses to well- 
established solutions (e.g. the diffusion limit) for end-member cases. 


Following Wordsworth and Pierrehumbert (2013), the diffusion of water through non- 
condensible background gases is given by 


Piip = uo 1,0 (1/H, = V9 ) 
by o. co, pCO, + by o-n, pN, + by oo, pO, 


Dio = 
pCO, + pN, + pO, 
8 3 14T s (26) 
LO = 
Hu,o8 
.314T. 
H, == 8 3 — strat 
HE 


Here, b, , (mol/m/s) is the binary diffusion coefficient of the i-th species through the j-th 


species (Marrero & Mason 1972; Zahnle & Kasting 1986). These are weighted by the 
stratospheric mixing ratios of each non-condensible constituent (CO2, N2, and O2), which 
are obtained from the atmospheric profile calculations (see above). The scale height of 
water, H„,o (m), and of the background gases, H, (m), depend on stratospheric 


temperature, 7;,,, = 200 K. The stratospheric water vapor mixing ratio is f, , and the 


strat 


diffusion limited escape flux is p, (mol HO0/m*/s). The diffusion-limited flux should 


arguably be set by the diffusion of atomic hydrogen through background gases since 
eddy diffusion dominates vertical transport at altitudes where molecular water is more 
abundant than atomic H and O (Catling & Kasting 2017), but our conservative approach 
minimizes oxygen accumulation from H escape. 


To calculate XUV-driven hydrodynamic escape of H, and associated O and CO» drag, we 
follow Odert et al. (2018) and Zahnle and Kasting (1986). The XUV-energy mass loss rate, 
Puy (kg/m^/s) is specified by the following equation: 
_ E(Fyyy X 96 >Etowxuv F xuv", 
XUV ^. 
4GM, 


o (27) 
Here, Fy, is the XUV flux (W/m?) received from the star (see stellar parameterization). 


The efficiency of hydrodynamic escape, e, is a function of atmospheric composition and 
XUV stellar flux, as described below. 


In general, the XUV-driven mass flux will be partitioned between H loss, O drag and, very 
under high XUV fluxes, CO» drag. The hydrogen escape flux, ®,, (molecules H/m*/s), can 
be obtained by analytically solving equations (4), (5), and (6) in Odert et al. (2018): 
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y moXo8(Mo -My )by_o | co, X co, 8 (Mco, — my )bu-co, 
a k,T(- Xp) kT 0+ by co Xo ! bo-co,) 


o 


Mco, X co, by co, gn, —-m,) Meo, X co, b, -CO, Xo! bo-co, 


(19 b, co, Xo / b, )k,T + Xo) It b, co, Xo / b; co, 


0-CO, 


(28) 
CO, 
Mco, X co, Meco, X CO, (b, -CO, Xo! bo-co, ) 


l+by -C0, Xo! bo-co, Ib, -CO, Xo! bo-co, 


-Oo, m, +MoXo + 


Here, m, (kg) is the mass of the i-th species, X, is the stratospheric mixing ratio of the i- 
th species, where we conservatively assume CO: is not dissociated to minimize the drag 
of carbon, T (K) is stratospheric temperature, k, is Boltzmann's constant, and 5, , is the 
binary diffusion coefficient of i through j. We refer the reader to the original paper for 
the details. Crucially, once the loss of hydrogen is known, the oxygen fractionation factor, 
Zo, can be obtained: 


zi gno -m,)b, , (29) 
o, k,T + X,) 


Xo 


If y, >0, then the hydrodynamic wind drags oxygen. The carbon dioxide fractionation 
factor can then be similarly calculated: 


LS 8(Mco, ~My )by co, /|(®,k,T) + by co, Xo d=%) + by co, XoXo I bo co, 
l+ by co, Xo ! bo, co, 
(30) 


Xco, = 


If %o, >9 then carbon dioxide is dragged along in the hydrodynamic wind and the 


corresponding escape fluxes (mol/m?/s) of O and CO» are given by the following 
expressions: 
0,-0,X,y, 


(31) 
Poo, = ® , X co, X co, 


Alternatively, if y, >0 and 7% o, <0, then oxygen is dragged but not carbon dioxide, and 
the corresponding escape fluxes are as follows: 
o, =0/(m,, +moXo%o) 
P, =(® -P my )/m, (32) 
Poo, =0.0 


If equation (29) yields y, <0, then then the XUV-driven escape flux is too small to drag 
oxygen. Carbon dioxide fractionation must be recalculated with y, = 0. If co, < 0.0 then 
the only escaping gas in hydrogen: 

o, = o/m, 


D, =P, =0.0 e 


CO; 
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In rare case where y, «0 and Yo, > 0, then loss rates as follows: 
Poo, =0,,%,X co, 
o, =0.0 


(34) 


It is convenient to convert these molecular escape rates to molar escape rates: 


Oxuv-y =P yMy [My 
9xuy-o = Pomo / uo (35) 


Pxuv-co, = Poo, Thco, / Hco, 


If the diffusion limited escape flux exceeds the XUV-driven H loss, then the diffusion 
limited flux is adjusted downwards: 


ug - min TON E /2 Pag) (36) 


Finally, we combine our expressions for diffusion-limit hydrogen escape and XUV-limited 
escape fluxes to obtain general expressions for the escape fluxes of hydrogen, Ey 


(mol/s), oxygen E, (mol/s), and carbon dioxide, Eco, (mol/s): 


W 20.9 +W. 
E, 1 12g 2xuv-u es 


w, +w, 
w. 
s, "(Suv lear? (37) 
w +w, 
W:Ọxuv-co, 2 
Eco, = Arr, 
w, +w, 


Here, the weightings w, and w, are a function of the H-abundances in the upper 
atmosphere: 
w =A (213-Xy) 


ks (38) 
w, = Xy 

These weighting functions ensure diffusion-limited H escape for low stratospheric 
abundances, and a smoothly transition to XUV-driven escape as the upper atmosphere 
becomes steam dominated. The precise transition abundance is unknown and will, in 
general, depend on conductive and radiative cooling of the upper atmosphere as well as 


downward diffusive transport. Here, it is represented by the free parameter, 4, which 


tra! 


ranges from 10? to 10° and is sampled uniformly in log space. 


The efficiency of hydrodynamic escape is parameterized by loosely following the 
approach of Wordsworth et al. (2018). If the XUV-stellar flux is insufficient to drag 
oxygen, then the efficiency is equal to a constant, é,,,,.,y, Which is randomly sampled 
from 196 to 3096. Alternatively, if the XUV-stellar flux exceeds what is required to drag O, 
then some portion of the excess energy, ¢ , goes into driving further escape, whereas 
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the rest, 1—¢ , is assumed to be efficiently radiated away. The efficiency factor, ¢ , is 


randomly sampled from 0-100% for complete generality. This leads to the following 
function for the efficiency of hydrodynamic escape: 


Enna Where E1wxuv xu" p " gm, — My b... 
4GM , k,T * X,) 
= b 
E(Fyuy > X 009s Etowxuv) = m 8 (Mo — My )by o 
i k,T(o* Xo) Elowxuv E xuv "p gno — My )by_o 

CE lowxuv + a =) RO VEM. a where ————— —— > Wy ————————— 

Fry, 4GM , koX 

4GM , 

(39) 


Fig. S6 illustrates our escape parameterization, with Monte Carlo outputs plotted for 
appropriate ranges in uncertain escape parameters. 


FXUV = 1 W/m2, H20/O2 atmosphere — fH2O = 90, fCO2 = 5, fO2 = 5 prcnt. 
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Fig. S6: Illustrative examples of the escape parameterization. Each line denotes a different 


calculation, where uncertain escape parameters ¢,,.,,,, and 4,,, have been randomly 


sampled. Black lines denote the cold trap diffusion limit (eq. (26)), red lines denote XUV- 
driven hydrodynamic escape (eq. (35)), and green lines show the weighted combination 
that is the escape parameterization in our model (eq. (37)). On the left-hand side, the 
stellar XUV flux is held constant, and escape fluxes are plotted as a function of the 
stratospheric H mixing ratio (for an atmosphere without CO»). On the right-hand side, 
the composition of the upper atmosphere is held constant, and escape fluxes are plotted 
as a function of the stellar XUV flux received by the planet. 


A.10) Solid-state evolution: Outgassing and crustal production 

The outgassing model follows that described in Wogan et al. (2020) where we calculate 
redox-dependent speciation of volatiles between melt and gas phase. Given mantle 
concentrations of H2O and CO: (by mass), corresponding melt fractions can be calculated 
assuming accumulated fractional melting: 


(i- a7)" ) m 


fh. HO = solid -H40 ái 
8 y M nante 
( = d = yn ) M ion 

ST neit-CO, = — ki 2 (41 ) 


Y M 


mantle 


Here y is the average melt fraction over the portion of the mantle where melting occurs 
(see below). These melt fractions, along with mantle oxygen fugacity (eq. (19)) and 
magma chamber pressure-temperature conditions are used as inputs to the outgassing 
calculations. For subaerial outgassing, the outgassing pressure is the pressure 
overburden of the atmospheric inventory, whereas for submarine outgassing, the 
pressure overburden is the atmospheric + ocean inventory. All outgassing is assumed to 
occur at the solidus temperature. Given these inputs, the outgassing thermochemical 
equilibrium model (Wogan et al. 2020) outputs gaseous mixing ratios, f,, for outgassed 
CO», H20, H2, CO, and CH,, as well as the moles of gas per total moles of melt plus 
gaseous species, œ. Note that the outgassing model does not consider the evolution of 


melt composition and oxygen fugacity along a degassing path; instead, we assume that 
the melt oxygen fugacity is buffered to that of the source rock, and that outgassed 
volatiles are determined by the equilibrium gas phase mixing ratios at surface pressure. 
Moreover, the molecular oxygen fraction of the degassed mixed is conservatively 
assumed to be negligible; the only possible source of atmospheric oxygen in the model 
is H escape. 


Gaseous mixing ratios can be convolved with melt production, MP , (m?/s, described 
below) to calculate outgassing fluxes, V, (mol/s), for each species: 
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fia Age. 


l-a, 


V,=MPxp, x ———- (42) 


Here, 9, = 15.5 mol magma/kg magma is the inverse molar mass of magma, which is 
assumed to be constant. The overall O2-consumption sink from outgassed volatiles is the 


summation of reducing species: 
V, uc = 0-5Vy, + 0.5Voo + Veen, (43) 


Overall outgassing fluxes are the combination of subaerial and submarine contributions, 
weighted by the surface land fraction, LF (see below): 
V y eden pe + V, Ptbmarine a —LF) (44) 


i 


To obtain mass fluxes, these molar outgassing fluxes must be weighted by their 
respective molecular masses: 


F, outgas-CO, — co, Veo, 


(45) 
tae = Hn,oVinno 
At any given time, the average melt fraction of freshly produced crust is given by 
integrating the melt fraction from the radius at which mantle temperature equals the 
solidus to the surface: 
y (T(r), Tie Eorsi T r), T sodes UE indo) 2 r)) x 4rr°dr 
y = r=r (T =T;otidus) — (46) 


A 
Azr^dr 


r=r(T=Tyotidus ) 


Here, the melt fraction at any given radius is given by the following expression: 


0, T(r) < T uds LC sete eg r) 
y Qn), T ola ? Ty idus) = EN l, T(r) > T iquidus( EA (47) 
T(r) = T aus (E easi ht 2 r) 


, otherwise 
T nd Caner ? r) = T olidus UP ririn ? r) 


Refer above for expressions for the solidus and liquidus. 


We assume the crust thickness, d.,,,, (m), is equal to the thickness of the melt layer, 


which can be calculated from the average melt fraction and total melt volume: 
1/3 


r=r, 


3 _ 
d. = f = r ——Ww Í Azr^dr (48) 
4 r-—r(T-T uius) 
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Once crustal depth is known, the melt production rate, MP (m?/s) can be calculated 
from the plate velocity, v,,,, (m/s), and ridge length, Liage (M): 


MP =l. _ xd. xv 


ridge crust plate 


plate 


(49) 


However, from the expression for half-space cooling of oceanic crust (e.g. Kite et al. 
2009; Turcotte & Schubert 2002, Sec 4.15), we can relate interior heatflow to ridge length 
and spreading rate as follows: 


2 
Shitty | Ey Fe Ye (50) 
2k (Ti x T, ) T 


Substituting this into eq. (49), melt production can be calculated from only the heatflow 
and crustal depth: 


2 
Arr,” 
MPs eee LES (51) 
2k (T. .-T,) Arr, 


sur) p 


Note that we are assuming the entire planetary area is involved in plate tectonics, which 
might cause us to underestimate melt production slightly. Plate velocity can be 
estimated by assuming a plausible ridge length, 3 zr, : 


eee MP | MP (52) 


plate Amr 


crust ridge crust p 


Plate velocity is only used for calculating seafloor weathering rates. Outgassing fluxes 
and crustal sinks of oxygen all depend on melt production. Note that for Venus and 
stagnant lid exoplanets the assumption of plate tectonics may overestimate melt 
production. Future versions of the model will explore a stagnant lid regime, but plate 
tectonics ought to maximize melt production and therefore maximize geological sinks of 
oxygen, which are the focus of this study. 


A.11) Solid-state evolution: Weathering 

Carbon is transferred from surface reservoirs to the interior via silicate weathering. 
Silicate weathering is the combination of continental weathering and seafloor 
weathering. To estimate this partitioning, we need to calculate the average depth of the 
ocean and corresponding land fraction. Following Cowan and Abbot (2014), we assume 
that there is a maximum ocean depth, d (m), above which there are no 


submerged continents. 
d s = a = Lipa rn )M pua - Ho 
d =11400(9.8/ g) 


ocean—max 


(53) 


Given these two quantities, we approximate the planetary hypsometric curve (proportion 
of land as a function of elevation) with a power law, and use this to calculate average 
land fraction, LF , and land fraction relative to the modern Earth, RLF: 
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LF =max {0.1 = (doean / TEA y^ 
LF ix 
1-(2.5/114) 


Clearly these are crude approximations, but they capture the fact the planets with a few 
Earth oceans ought to have some subaerial land, but that for large water inventories all 
crust is submerged. In any case, total land fraction does not have a large impact on 
weathering feedbacks (Abbot et al. 2012). 


Continental weathering fluxes, W 


cont 


(kg/s), and seafloor weathering fluxes, W,,, (kg/s), are 
given by expressions similar to those described in Krissansen-Totton et al. (2018): 


Wyr — Woef E 10 9308-7) V plate exp — E, 1 RS 1 (55) 
4 3cm/ yr 8314| Tep 285 
7 T „—285 
Wont = Woef x RLF x pCO, exp a=] (56) 
; 350 ppm Toa 


Here, the multiplicative factor W. 


soet = 4000 kg/s is chosen to approximately reproduce 
modern Earth fluxes (or rather is the weathering flux required to balance mantle-derived 
CO» outgassing). Unknown, randomly sampled parameters include the temperature 


dependence of continental weathering, Tpu = 5-30 K and the CO? dependence of 


continental weathering, y = 0.1-0.5. The temperature dependence of seafloor 
weathering, E, = 90 kJ/mol, and pH-dependence of seafloor weathering are assumed to 
be known constants. More sophisticated weathering parameterizations that account for 
kinetic dependencies, thermodynamic solute concentration limits, and a precipitation- 
limited runoff dependence have been proposed (Graham & Pierrehumbert 2020), but 
given the uncertainties in geological parameters that feed into such models, our simple 
CO» and temperature dependent formulation is adequate for providing a crude 
thermostat. 


Total weathering will be the summation of continental and seafloor weathering, 
weighted by the fraction of liquid water at the surface, (1— fr, ,, ,) . This ensures that 


tmo-H, 
weathering tends to zero as oceans evaporate. We also include a possible supply limit to 
weathering as an unknown variable that could potentially limit the rate at which dense 
pCO: atmospheres are sequestered if the supply of erodible rock is low. This supply limit, 
W, p-m = 10° to 10" (kg/s), is not coupled to crustal production since not all newly 
produced crust will necessarily be delivered to the surface to be weathered. The overall 
expression for CO2-removal via weathering is therefore: 


Fveather-COy -min (W. i , Jiquid-H,0- fraction (Woon + Wr )} (57) 
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Because we are modeling a plate tectonics regime, we assume that all carbon dioxide 
removed by weathering is returned to the mantle. This can be seen in equation (16), 
whereby the weathering flux, Fyeaner—co,, transfers carbon from the surface fluid reservoir 


to the interior mantle reservoir. 


A.12) Solid-state evolution: Deep water cycle 


The transfer of water from the mantle to the surface is already specified by F, In- 


utgas-H,0 * 
gassing of water from the surface to the interior is controlled by water-rock reactions 
and tectonic processes. There is considerable uncertainty in the function dependencies 
of ingassing reactions (e.g. Cowan & Abbot 2014; Komacek & Abbot 2016; Schaefer & 
Sasselov 2015). Here, we loosely follow the approach in these papers, but introduce 
several unknown parameters to account for different deep hydrological cycle feedbacks. 


The maximum depth in the crust within which hydration reactions may occur is 
approximated by: 


DT fomes gc. 
dos E k p surf surf = si surf ) (58) 
T,-T, Im 


In surf 
Here, 973 K is the maximum surface temperature for serpentine stability (Schaefer & 
Sasselov 2015), and k = 4.2 W/m/K is the thermal conductivity of silicates (c.f. equation 
(4)). Since we assume no hydration occurs below the crust, it is also helpful to define the 
fractional depth of hydration as the ratio of the hydration depth to the crustal depth, or 
1.0 (whichever is smaller): 


Írydr-depin = min ae [ders > 1.0} (59) 


Water loss from surface reservoirs can be conceptually partitioned into hydration 
reactions that add water to the solid interior but do not alter atmospheric redox state, 
e.g. 

2(Fe,Mg,Ca)O+H,O — (Fe,Mg,Ca)(OH), (60) 


and hydration reactions that oxidize the solid interior and result in outgassed hydrogen 
(i.e. hydrogen that is lost to space under anoxic conditions, or recombines with 
atmospheric oxygen under oxidizing conditions): 

2FeO+H,O > H,+2FeO, , (61) 


Hydration reactions that add water to the interior, equation (60), are parameterized as 
follows: 


d d M solid -H,0 
Fy o iyd: = min = ,l.0 Snydr—depth Pa frae MP Pm max 0,1 = : 


max-—ocean solid -H,O—max 


(62) 


Here, Jm rur = 10? to 0.03 (sampled uniformly in log space) is the unknown efficiency 


frac 


of hydration reactions. We are assuming that, at most, hydrated crust is 396 water by 
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mass (Schaefer & Sasselov 2015), but this number could be much less depending on the 
tectonic regime and cracking of the crust. Some unknown portion of this crustal water 
may be returned to the surface via arc volcanism, for instance. After considering this 
efficiency factor, we assumed all hydrated crust is returned to the mantle. We assume a 
linear dependence on ocean depth for as long as there is emerged land, then no ocean 
depth dependence beyond this. Additionally, we assume that the return of water to the 
interior tapers off as the water content of the mantle approaches its maximum value, 

M i 4,9. (K9). This unknown variable is sampled randomly from 0.5-15 Earth oceans 


(Komacek & Abbot 2016). 


Hydration reactions that oxidize the crust and remove water from surface, but do not 
add water to the interior (i.e. serpentinizing reactions that produce hydrogen, equation 
(61)) are parameterized as follows: 


d M. 
Fy,o-ser = fu —oxid min [ema „MP p, Xy, So et 
P i max—ocean i i M sia ro + M uia reo, 
(63) 
Here, faoa is another unknown efficiency parameter (10? to 10°') representing the 


fraction of crustal iron that is oxidized via hydration reactions (Lécuyer & Ricard 1999). 
This efficiency parameter represents the degree of serpentinization in water-rock 
reactions. Total ingassing contributions are calculated as follows: 

=F, +F, eee 


ingas—H,O-loss H,O-hydr H,O-serp 
l 29 j ST n (64) 


Edad ud Fy o hyar 


The former, E 


ingas—H,O-loss ! 


represents the flux of water lost from the surface reservoir, 
whereas the latter, F,,,. |, ,,;,, represents the flux of water into the interior via 
subduction of hydrated crust. This transfer of water can also be seen in equation (16). 
Note that these loss and gain quantities are not necessarily equal in because, as 


discussed above, H2 produced by serpentinization may be lost to space, thereby 
permanently removing water from the planet. 


A.13) Solid-state evolution: Planetary redox budget 

The interior may become oxidized via outgassed of reduced species (discussed above), 
dry oxidation, and wet oxidation. First, we consider dry (direct) crustal oxidation, which 
can be represented by the following reaction: 


2FeO+0.50 , — 2FeO, , (65) 
The flux of this dry crustal sink is parameterized by the following equation: 
M u 
Fes -osid pz (1 ix RLF) LN X ke solid Feo DMP s oe (66) 
: : M via rio i M oia reo, 2 FeO 
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Here, there is a land fraction dependence to ensure that no dry oxidation of the crust 
occurs when the surface is covered in water. The unknown efficiency parameter f, xia 


(10% to 10%) is discussed in the main text, and is the fraction of reduced iron in newly 
produced crust that is oxidized. 


We also modify the melt production term to represent the fact that there is some 
maximum amount of surface melt accessible to oxidation via diffusion through extrusive 
lava flows. The diffusivity of oxygen in silicate melts is D. « 10°’ cm?/s (Canil & 
Muehlenbachs 1990), which implies a downward diffusion depth of ~./D,,r per unit 
time, t. The maximum downward oxygen flux averaged over one year is therefore given 
by: 


, XD, x3.15x10's 

le 

P 3.15x10’sx100cm/m (67) 
= Sizi x 2.87 x10° m /s = dus x9.05 x 10° km? / yr 


E có; = Joni 4r 


Here, the efficiency factor, f, is the average fraction of the planetary surface that is 
continuously molten due to extrusive volcanism after the magma ocean has solidified. 
Because thermal diffusivity exceeds chemical diffusivity—and because mean surface 
temperature is below the solidus by definition—extrusive magmas will form a low 
permeability solid crust as they cool, precluding continuous diffusion of oxygen. Thus, 
even for extreme rates of resurfacing due to high heatflow, fa is likely low. For 
example, on lo, where average internal heatflow is 1-3 W/m? (Veeder et al. 2012), only a 
few km? of the surface is estimated to be molten at any given moment (Mura et al. 2020). 
We uniformly sample f,,,, = 10 to 1 in log space for full generality. Given this upper 
limit on melt oxidation, the amount of oxidizable melt is given by: 


$a DER 


domat 


MP, 


surface = min L (68) 
This avoids unrealistic oxygen sinks during the transition from magma mush to solid 
mantle where melt volumes are extremely high, but the melt accessible to atmospheric 
oxygen via diffusion through extrusive magmas limited. The value for f£, does not 
affect any of our model scenarios except for Scenario 3, where values exceeding 0.01 are 


required for oxygen accumulation (Fig. S1d). 


The wet oxidation of the interior from hydration reactions is already the appropriately 
weighted serpentinization flux: 


Fonts rs Fa,0-serp (69) 
FeO 


Recall the term Fia- ,,,, 


atmosphere-ocean reservoir to the interior, whereas the term F,,;, quia represents the 


represents the total flux of free oxygen lost from the 
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total flux of free oxygen gained by the interior. Oxidized crust is assumed to be mixed 
back into the mantle on long timescales (equation (16)) via subduction, or slab 
delamination. We also assume that the outgassing of reduced species must ultimately be 
balanced by oxidation of the crust, and so the so net oxidation of the interior equals the 
reduction of the fluid reservoir: 


| ac fuia a IT osi * PE eod + Ho, Vo, -sink T P dd (70) 


The only exception to this equality is when atmospheric oxygen levels are very low, and 
fluxes are modified for numerical reasons (see below). 


Text B. Numerical approach 

All code is written in python and available open source ([DOI TBD]). The system of 
differential equations was solved explicitly with either RK45 or RK23 using the solve ivp 
module in scipy. The maximum timestep was shorter for the magma ocean phase (10000 
years) compared to the temperate evolution phase (10° years). 


To avoid sawtoothing and excessive computation time at low reservoir abundances, 
various adjustments were made to the differential equations to ensure adequate 
performance. First, if atmospheric CO» dropped below 50 Pa and if weathering exceeds 
outgassing, then the time-derivative of surface carbon dioxide was set equal to zero. This 
may mean climate evolution is slightly inaccurate at low pCO», but the effects are minor. 
Second, atmospheric oxygen is similarly prevented from dropping below 0.1 Pa. If 
atmospheric oxygen is below 0.1 Pa and if oxygen production via escape is less than 
oxygen consumption, then the following adjustments are made. The atmosphere is 
assumed to be in an anoxic steady state and so oxygen atmospheric sinks are set equal 
to the escape source. However, since oxygen production via escape is less than oxygen 
sinks (outgassing and Hz, wet and dry oxidation reactions), then the interior is assumed 
to be oxidized by the difference as the excess reductants (H2) escape to space. This might 
slightly overestimate mantle oxidation because some reductants (e.g. CO) will get 
photochemically oxidized rather than escape, but the redox budget of the atmosphere is 
not directly affected, and the effects on mantle redox evolution are negligible. 


Text C. Venus Model validation: 

To further validate the model, we demonstrate that it can broadly reproduce the known 
atmospheric evolution of Venus. To model Venus, all parameters are kept the same as 
for Earth except for planet radius, mass, planet-star separation, and albedo 
parameterization (see Table S2). Fig. S7 shows all model outputs that reproduce modern 
Venus, which is defined to be atmospheric oxygen < 107? kg («0.2 mbar), atmospheric 
CO» exceeding 2x10? kg (40 bar), no surface water ocean, and atmospheric H2O < 
2x10'6 kg (~3 mbar). 
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S7: Model runs that reproduce modern Venus conditions. Note that there are two 


qualitatively different histories that can reproduce modern Venus (top left). Either Venus 
was always in a runway greenhouse phase and never condensed a surface ocean, or 
Venus maintained a temperate surface for several Gyr before transitioning back to 


runaway greenhouse as solar insolation increased. 


Note that our Venus model somewhat overpredicts modern day Venusian heat flow and 
melt production because we assume a plate tectonics model (Nimmo & McKenzie 1998). 
We save a more detailed comparative study of solar system planets with stagnant lid 


tectonics and resurfacing events for future study. 


In Fig. S8 we plot key parameter values for model runs that 


successfully reproduce 


modern Venus. Initial volatile inventories are likely small (Fig. S8a) and dry crustal 


oxidation must be relatively efficient (Fig. S8c). 
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Fig. S8: Parameter values for individual model runs in Fig. S7. These shows the require 
initial volatile inventories (bottom left) and dry crustal oxidation efficiency (top left) 
required to reproduce modern Venus. A broad range of XUV escape efficiencies (top 
right) and fractional molten areas (bottom right) are permissible. 


Text D. Hydrous Mantle Sensitivity Test 

Here, we test whether modifying the solidus for hydrated mantles affects oxygen 
accumulation in the waterworlds scenario. Following Katz et al. (2003) we modify our 
expressions for the solidus and liquidus as follows: 


_T(r)exp(10°(-r, +r+ 10°)) T, G)exp(10? (r, -r -10°)) 


T sius T> Pover urden^ ~~ 
E exp(10°(-r, +r +10°)) +exp(10°(r, -r —10°)) 


overburden 


T, = max fi 170, 1420 +104.42x10° gp, (r, -r)+104.42x10° P = 4.74108 (M uu io [M manie yt 


overburden 


T, = max 11170, 1825 + 26.53x10? go. (r —r)+26.53x10°P —434x10* (M... /M.) - 
2 § Pn p solid CHO mantle 


Ty idus (r, Poverburaen) = T otidus (r, Tousiasdós ) t 600 
(71) 


This modification is a crude polynomial fit to Fig. 3 in Katz et al. (2003). It accounts for 
the depression of the solidus as the water content of the mantle increases, but the 
solidus never drops below about 900°C, which is the approximate water saturation limit. 
Waterworld calculations repeated with this hydrated solidus are shown in Fig. S9. 
Although crustal production is elevated compared to the anhydrous solidus case in the 
main text, outcomes are qualitatively similar to Fig. 5 in the main text. Significant abiotic 
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oxygen accumulation becomes increasingly likely for initial water inventories exceeding 
100 Earth oceans (Fig. S10). 
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Fig. S9: Same as Fig. 5 in the main text except the solidus decreases with mantle 
hydration. Early crustal production is elevated (bottom left), but outcomes are 
qualitatively similar. Oxygen sinks are shut down by the pressure overburden after a few 
billion years and oxygen accumulates. 
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Fig. $10: Same as Fig. 6 in the main text except the solidus decreases with mantle 
hydration. Abiotic oxygen accumulation is somewhat less frequent, and occurs at higher 
initial water inventories, but results are qualitatively the same. 
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Text E. Impact Ejecta O2-sinks 


Here, we test whether the delivery of reducing materials from impactors could 


potentially draw down oxygen produced in desertworld scenarios. We introduce an 


impactor flux, F,,, (kg/yr), that diminishes exponentially with time: 
PB = E exp(t/t,,,,, ) 


Here, the coefficient F° 


imp 


(72) 


is randomly sampled (in log space) from 10"' to 10? kg/yr and 


the decay time, taa (Gyr), is randomly sampled from 0.06 to 0.14 Gyr. These ranges are 


adopted because they approximately reproduce plausible estimates for impactor fluxes 
in the Hadean and early Archean, both with and without a late heavy bombardment 
(Kadoya et al. 2020). Additionally, we assume that impactors are 3096 metallic iron by 
mass, and that 10096 of this iron is completely oxidized to ferric iron instantaneously 
upon impact, depleting atmospheric oxygen. Fig. S11 shows our desertworld calculations 
repeated with this impactor flux. We find oxygen accumulation and retention for several 
Gyr is still possible, although only when the total impactor flux is low. Fig. $12 shows 
oxygen accumulation after 4.5 Gyr as a function of total impactor flux. Abiotic oxygen 


may accumulation for impactor mass fluxes < 10% kg. 
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Fig. S11: Same as Fig. 7 in the main text except impact ejecta sinks for oxygen have been 


added. Retention of abiotic oxygen is still possible if impactor fluxes are low. 
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Fig. S12: Scenario 3 abiotic oxygen accumulation as a function of total impactor flux. 
Large impactor fluxes preclude the retention of abiotic oxidation if all impactor material 
is efficiently oxidized. 
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Text F. Stratospheric Temperature Sensitivity Test 
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Fig. S13: Sensitivity ot talse positive results to stratospheric temperature. 
Calculations in main text were repeated, but stratospheric temperature is a free variable 
that is randomly sampled from 150 K to 250 K. Final oxygen accumulation after 4.5 Gyrs 
is plotted as a function of stratospheric temperature. Subplot (a) shows all Scenario 1 
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false positives model runs, whereas (b) shows all Scenario 2 false positives, and (c) shows 
all Scenario 3 false positives. In (b) red dots denote leftover oxygen after the magma 
ocean, whereas blue dots show model runs where all oxygen produced during the 
magma ocean phase is sequestered in the mantle, and oxygen builds up subsequently, 
as described in the main text. 


Fig. $13 shows the sensitivity of each false positive scenario to stratospheric temperature. 
Each dot is a model run representing an oxygen false positive. For the first scenario (Fig. 
$13a), abiotic oxygen only occurs when stratospheric temperature exceeds ~200 K. This 
is because, at lower temperatures, the cold trap becomes more effective and H escape 
(and therefore O accumulation) is throttled. 


For the waterworld scenario (Fig. S13b) oxygen accumulation may occur at any 
stratospheric temperature. However, this is more probable—and abiotic oxygen 
abundances are greater—at higher stratospheric temperatures. On waterworlds, cold 
stratospheres are not necessarily expected because an N2-dominated atmosphere with 
low CO» is a probable outcome (Fig. 5), especially if continuous CO2-drawdown via 
weathering occurs (Nakayama et al. 2019). Moreover, modest oxygen accumulation 
would result in ozone formation, that would further warm the stratosphere, potentially 
resulting in a positive Oz-accumulation feedback that is not considered here. Note that 
there are two subclasses of oxygen false positives in Fig. S13b, denoted by red and blue 
dots. The blue dots show model runs where oxygen accumulated during the initial 
magma ocean is completely sequestered in the mantle upon magma ocean solidification, 
whereas red dots denote scenarios whereby appreciable oxygen is left over after magma 
ocean solidification due to the high pressure-temperature conditions of the overburden 
suppressed solidus. 


The third, desertworld scenario (Fig. S13c) is largely insensitive to stratospheric 
temperature. This is because water loss and oxygen accumulation occur immediately 
after magma ocean solidification while the steam-dominated atmosphere persists. There 
is no effective cold trap in the steam-dominated atmosphere and so escape fluxes are 
insensitive to stratospheric temperature. 


Text G. Reducing Mantle Sensitivity Test 


The nominal model assumes Earth-sized planets undergo rapid core formation with 
mantles that quickly approach ~FMQ+2. While this is a common assumption when 
modeling magma ocean evolution of the early Earth (e.g. Zahnle et al. 2007; Zahnle et al. 
2010) and of terrestrial exoplanets (e.g. Schaefer et al. 2016), this may not be true for all 
terrestrial planets. To investigate the effects of a more reduced initial mantle, sensitivity 
tests were performed whereby the initial magma ocean was endowed with a smaller 
amount of free O (0.5x 10?! to 1.5x 10?! kg), such that after 4.5 Gyrs of evolution, mantle 
oxygen fugacity is closer to the iron-wüstite buffer than FMQ. Additionally, following 
Ortenzi et al. (2020), we modified the outgassing model such that the melt-solid 
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partitioning of carbon is controlled by graphite saturation. The concentration of carbon 
dissolved in a graphite-saturated melt in redox state dependent: 


X co, ,grahite—sat = [E qu NTE N) / (1 * Rian oada) 


(73) 
T (44 / 36.594) X cO. gihe sii / (1 z (1 z 44 / 36.594) X co, grahite-sat ) 


CO, , grahite—sat 


Here, fO, is mantle oxygen fugacity, and we are converting between dissolved 


carbonate and carbon dioxide concentrations. The temperature and pressure-dependent 
equilibrium constants are defined as follows: 
log, (K ) = 40.07639 — 2.53932 x10 7 T + 5.27096 x 10 * T^ 0.0267 (P-1 )/ T 


log, (K -6.24763 — 282.56 / T — 0.119242 (P —1000)/ T 
(74) 


1, graphite 


annie) = 


To calculate melt concentrations for outgassing calculations, we take the minimum of the 
concentrations in equations (41) and (73): 


par Vkeo, 
(1 = ad d | M ii co, (75) 
W M ade 


foco; = min X co, grahite—sat ? 


Taking the minimum ensures that graphite saturation does not overestimate dissolved 
carbon concentrations under oxidizing conditions and when the total carbon content in 
the mantle is low. 


Finally, while we do not explicitly account for graphite precipitation during magma ocean 
solidification, we set kco, =1.0 in equation (16) to allow for greater retention of carbon in 


the mantle. It should be emphasized that these modifications do not constitute a fully 
consistent model of reduced mantle planetary evolution because the radiative transfer 
model does not allow for CO and H2 dominated atmospheres. However, for post magma 
ocean evolution they are adequate approximations. 


Fig. S14 is identical to Fig. 2 in the main text except for the reducing mantle initial 
conditions and other assumptions described above. Once again, an anoxic atmosphere is 
assured after 4.5 Gyrs of evolution because crustal oxygen sinks overwhelm oxygen 
sources. Fig. S15 is the reduced mantle equivalent of Fig. 3 in the main text showing high 
CO2:H20 oxygen false positives. This scenario is largely unchanged by a reducing mantle; 
magmatic outgassing does not occur due to the high pressures and low mantle volatile 
concentrations following magma ocean solidification. Gradual oxygen accumulation may 
occur after several Gyrs of H loss to space. Fig. $16 is the reduced mantle equivalent of 
Fig. 5 in the main text showing waterworld oxygen false positives. The pressure 
overburden of a large surface ocean once again shuts down crustal production after ~1 
Gyr, thereby removing all crustal oxygen sinks and permitting atmospheric oxygen to 
accumulate. Fig. S17 is the reduced mantle equivalent of Fig. 7 in the main text showing 
desertworld oxygen false positives. Although Scenario 3 is apparently unchanged by 
having a lower mantle redox, a fully self-consistent model that accounted for the high 
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CO-H?» content of the originally degassed atmosphere would likely preclude early O2 
accumulation, in practice. 
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Fig. S14: Nominal Earth evolution with a more reduced initial mantle. This is identical to 
Fig. 2 in the main text except (i) the initial free oxygen of the mantle is lower (0.5- 

1.5x 10?' kg), (ii) graphite saturation in reduced melts is accounted for, (iii) and carbon is 
partitioned into the solid phase during magma ocean solidification. 
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Fig. $15: Scenario 1 oxygen false positives with a more reduced initial mantle. This is 
identical to Fig. 3 in the main text except (i) the initial free oxygen of the mantle is lower 
(0.5-1.5 x 10?! kg), (ii) graphite saturation in reduced melts is accounted for, (iii) and 
carbon is partitioned into the solid phase during magma ocean solidification. The 
terminal magma ocean becomes more oxidized than the solid mantle as H escape 
oxidizes the combined melt-volatile reservoir faster than solidification transfers oxidized 
material to the mantle. Oxygen accumulation may occur after several Gyr because 
outgassing of C-bearing volatiles is negligible from the graphite-saturated mantle. 
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Fig. S16: Scenario 2 oxygen false positives with a more reduced initial mantle. This is 

identical to Fig. 5 in the main text except (i) the initial free oxygen of the mantle is lower 

(0.5-1.5 x 10?! kg), (ii) graphite saturation in reduced melts is accounted for, (iii) and 

carbon is partitioned into the solid phase during magma ocean solidification. Oxygen 


sinks are suppressed by the large pressure overburden of the surface ocean, the same as 
in the nominal oxidized-mantle calculations. 
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Fig. S17: Scenario 3 oxygen false positives with a more reduced initial mantle. This is 
identical to Fig. 7 in the main text except (i) the initial free oxygen of the mantle is lower 
(0.5-1.5 x 10?! kg), (ii) graphite saturation in reduced melts is accounted for, (iii) and 
carbon is partitioned into the solid phase during magma ocean solidification. The 
terminal magma ocean becomes more oxidized than the solid mantle as H escape 
oxidizes the combined melt-volatile reservoir faster than solidification transfers oxidized 
material to the mantle. While the persistence of oxygen is permitted in these 
calculations, in practice, the early atmosphere is likely too reducing to permit such 
oxygen accumulation. 


Text H. Stellar Separation Sensitivity Test 


While this study is not an exhaustive exploration of the oxygen false positive parameter 
space, nominal model calculations were repeated at 1.3 AU to show that oxygen 
accumulation is not dependent on being close to the inner edge of the habitable zone. 
Fig. $18 shows all Scenario 2 false positives at 1.3 AU. The increased stellar separation 
results in lower H escape fluxes, but oxygen accumulation may still occur if crustal sinks 
are suppressed by pressure overburden. Similarly, Fig. S19 shows all Scenario 3 false 
positives at 1.3 AU. Oxygen accumulation on desertworlds can similarly occur at larger 
stellar separations because oxygen accumulation occurs early during the steam- 
dominated atmosphere. Scenario 1 false positives do not occur at large stellar 
separations because even under a high CO? atmosphere, the runaway greenhouse state 
is not maintained after magma ocean solidification. 
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Fig. S18: Identical to Fig. 5 in the main text except the assumed planet-star separation is 
1.3 AU. Abiotic oxygen accumulation is still permitted due to overburden pressure 
suppressing oxygen sinks. 
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Fig. S19: Identical to Fig. 7 in the main text except the assumed planet-star separation is 
1.3 AU. Abiotic oxygen accumulation is still permitted due to early oxygen accumulation 
during the steam-atmosphere phase. 
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Text I. The Effect of Sulfur Outgassing 


Sulfur outgassing and burial may have played an important role in the oxygenation of 
Earth's atmosphere (e.g. Gaillard et al. 2011; Olson et al. 2019). While including a 
complete model of sulfur cycling is beyond the scope of this study, we present 
calculations showing that Earth-like sulfur mantle abundances are unlikely to 
qualitatively change our conclusions. 


Following Gaillard and Scaillet (2014), sulfur speciation is added to our outgassing model 
by adding the following system of equations to those already described in Wogan et al. 
(2020): 

In(x, ) 2 0.5x In(p; ) -0.5x1n( fo, ) -15.274 


Xs, -rota P Ex (Ps, + 2 Pus + 2Pso, gas + (1 Ag )xs,P 
In(Pso,) = In K; fo, )+ 0.5 In(ps, ) 
0.5 x In(p, ) =—In(Py,o) + W(K; fo, ^) + (Pu, 5) 


(76) 


Here, p, refers to the gas phase partial pressure of species n, x, _,,,., is the total sulfur 

mass fraction of the original melt, a, is the gas phase mole fraction, fọ, is the oxygen 

fugacity of the melt. The equilibrium constants K, and K, are defined as follows: 
K, = exp(1.90560415x10* / T —0.860366131) 


(77) 
K, = exp(4.35250424 x10* / T — 8.80403494) 


The equations in (76) are solved simultaneously with the system of equations describing 
carbon, oxygen, and hydrogen gas-melt speciation (Wogan et al. 2020) to determine 
total outgassing fluxes, and the overall outgassing oxygen sink: 

Vo, -sine 7 0.54, +0. Veo + Vey, € 15V, + V5, (78) 
Rather than attempt to model complete sulfur cycling, model outputs from Fig. 2 were 
post-processed to recalculate outgassing fluxes and oxygen sinks using the above 
model, assuming constant 300 ppm mantle sulfur abundances (von Gehlen 1992). Fig. 
S20 shows the results from this calculation comparing oxygen sinks with and without 
sulfur outgassing. The impact on total oxygen sink fluxes is relatively minor. Fig. S21 
shows the same calculations performed on the Scenario 1 false positive outputs from Fig. 
3. In this case, the high pressure for the dense CO» atmosphere prevents the exsolution 
of sulfur-bearing gases from the melt, and so total oxygen sink fluxes with and without 
sulfur are barely distinguishable. Moreover, for Scenario 2, the shutdown of crustal 
production inhibits all crustal sinks, regardless of volatile content. 
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Fig. $20: Comparison of oxygen sinks with and without sulfur-bearing outgassed 
volatiles. Subplot (a) shows outgassing sinks in the nominal model from Fig. 2 (red) 
compared to outgassing sinks if sulfur-bearing species are included (blue) and a 
constant 300 ppm mantle sulfur concentration is assumed (von Gehlen 1992). Subplot (b) 
shows total oxygen sinks with (blue) and without (red) sulfur-bearing volatiles. While the 
inclusion of sulfur-bearing species may result in a slightly larger oxygen sink, for Earth- 
like mantle concentrations the effect of sulfur outgassing is minimal. 
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Fig. S21: Comparison of oxygen sinks with and without sulfur-bearing outgassed 
volatiles. Subplot (a) shows outgassing sinks for Scenario 1 false positives from Fig. 3 
(red) compared to outgassing sinks if sulfur-bearing species are included (blue) and a 
constant 300 ppm mantle sulfur concentration is assumed (von Gehlen 1992). Subplot (b) 
shows total oxygen sinks with (blue) and without (red) sulfur-bearing volatiles. Including 
sulfur-bearing species has a negligible effect on outgassing sinks because the pressure 
of the dense CO» atmosphere inhibits exsolution of sulfur gases. 
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Table $1. Monte Carlo analysis and uncertain parameter ranges (for nominal Earth 


model). 
Nominal References/Notes 
range 
Initial Water 10?'-10? kg* | Approximately 1 — 10 Earth 
conditions oceans 
Carbon dioxide 107-107 kg* | Approximately 20-2000 bar (if 
no other atmospheric 
constituents). 
Radiogenic inventory (relative 0.33-3.0 Scalar multiplication of 
Earth) inventories in Lebrun et al. 
(2013) 
Mantle free oxygen 2x10*!- This ensures a post- 
6x10% (kg) solidification mantle redox 
around Quartz-Fayalite- 
Magnetite buffer. 
Solar Early sun rotation rate (relative | 1.8-45 (Tu et al. 2015) Implies XUV 
evolution modern) saturation time of 6 — 226 Myrs. 
and escape | Escape efficiency at low XUV 0.01-0.3 See escape section. 
parameters | flux, e. 
Transition parameter for cold- | 10°? — 10? See escape section. 
trap diffusion limited to XUV- 
limited escape, A, 
XUV energy that contributes to | 0-100% See escape section. 
XUV escape above 
hydrodynamic threshold, ¢ 
Carbon Temperature-dependence of 5-30 K (Krissansen-Totton et al. 2018) 
cycle continental weathering, Tpu 
parameters | CQ,-dependence of continental | 0.1-0.5 (Krissansen-Totton et al. 2018) 
weathering, » 
Weathering supply limit, W,,,, jim | 10°- 10" kg/s | (Foley 2015) 
Ocean calcium concentration, 10^-3x10" | See text for explanation (Halevy 
[d mol/kg & Bachan 2017; Kite & Ford 
2018) 
Ocean carbonate saturation, Q | 1-10 (Zeebe & Westbroek 2003) 
Interior Mantle viscosity coefficient, 10'-10?Pas | Fit to modern heatflow and melt 
evolution Vos production (see Fig. S3) 
parameter 
Crustal Crustal hydration efficiency, 10? to 0.03 Upper limit wt ?6 H2O in oceanic 
sinks Jac ne crust. Lower limit hydration 
oxygen and limited by cracking. 
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hydrological | Dry oxidation efficiency, f,.,.@ | 10% to 10% Plausible range of processes for 
cycle l Venus (Gillmann et al. 2009) 
parameters | Wet oxidation efficiency, 10? to 107 Based on oxidation of Earth's 
pnr oceanic crust (Lécuyer & Ricard 
1999). 
Maximum fractional molten 10^ to 1.0 Refer Supplementary Text A.13 
area, fro, 
Max mantle water content, 0.5-15 Earth (Cowan & Abbot 2014) 
M ta n0 max Oceans 
Albedo Hot state albedo, A, 0-0.3 See albedo parameterization. 
parameters | Cold state albedo, A. 0.25-0.35 See albedo parameterization. 


* We apply the additional constraint that the initial water inventory is greater than the 
initial carbon dioxide inventory, in accordance with typical carbonaceous chondrite 
abundances. 


Table S2: Changes in Monte Carlo parameters for different abiotic oxygen scenarios and 
Venus validation: 


Nominal High Waterworld | Desertworlds | Venus 
range COz:H20 S 
Initial H20 1071-10” kg* | 10*°-10°4 107-1075 10193-10707 10°°-10* kg 
cond. kg** kg kg 
CO> 1029-1022 kg* | 1 020-102? 1029-10225 10195-10205 1029-1022 kg 
kg** kg kg 
Albedo | Hot, A,, | 0.0-0.3 0.0-0.3 0.0-0.3 0.0-0.3 0-0.3 
param. | Cold, 0.25-0.35 0.25-0.35 0.25-0.35 0.25-0.35 0.2-0.7 
Ac 


* We apply the additional constraint that the initial water inventory is greater than the 
initial carbon dioxide inventory, in accordance with typical carbonaceous chondrite 
abundances. 

** High COz:H;O runs are a subset of this range with initial CO2/H2O>1 (by mass). 


Table S3: All fixed parameters used in the model. 


Parameter Value 

Planetary iron content (silicate mole fraction), 0.06 

Xe 

Average silicate density, p,, 4000 kg/m? 

Planetary radius, P 6371 km (Earth) or 6052 km (Venus) 

Core radius, r, 3460 km (Earth) or 3230 km (Venus) 

Planet mass, M, 5.972 x 10 kg (Earth) or 4.867x 107^ kg 
(Venus) 

Latent heat of fusion of silicates, H „sion 4x10? J/kg 
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Specific heat of silicates, c, 


1200 J/kg/K 


Thermal expansion coefficient for silicates, c 2x10° K1 
Critical Rayleigh number, Ra, 1100 
Thermal conductivity of silicates, k 4.2 W/m/K 
Thermal diffusivity of silicates, x 10° m?/s 
Convective heatflow exponent, 5 1/3 


Molecular mass of i-th species, x, 


Various (kg/mol) 


Crystal-melt partition coefficient for water, 
k 


H,0 


0.01 


Crystal-melt partition coefficient for carbon 
dioxide, kco, 


2x10? 


Binary diffusion coefficient of the i-th species 
through the j-th species, b, , 


Various (mol/m/s) 


Stratospheric temperature, T 200K 
Inverse molar mass of magma, 39, 15.5 mol/kg 
Weathering multiplicative coefficient, W, 4000 kg/s 
Activation energy seafloor weathering, E,,. 90 kJ/mol 


Mass of i-th species, m, 


Various (kg) 


Planet-star separation, D 


planet—star 


1.49610" m (Earth), 1.047x10!! m (Venus) 
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